3.7 Équation de diffusion non linéaire

3.7.1 Problème physique: diffusion dans une flamme

On veut calculer la répartition de température transversale dans une flamme (figure 3.32). Pour cela on utilise un modèle simple, en supposant que dans la direction transverse (notée $x$ ) on a des échanges de chaleur uniquement par diffusion et par rayonnement.

Figure 3.32: Profil de température transversal dans une flamme

Avec ces hypothèses, l'équation de conservation de l'énergie à l'état stationnaire dans la direction $x$ s'écrit:

$$\underbrace{-\frac{\partial}{\partial}\left(\lambda\frac{\partial... ...ce{\sigma(T^{4}-T_{0}^{4})}_{\mbox{rayonnement}}=\underbrace{Q}_{\mbox{source}}$$ (3.66)

Dans ce modèle, l'énergie $Q$ produite dans la flamme de largeur $2\delta$ par réaction chimique est diffusée par conduction et rayonné vers l'air extérieur à température $T_{0}$ . Pour le rayonnement, nous avons adopté un modèle simple de rayonnement de corps noir proportionnel à $T^{4}$ avec une constante de rayonnement $\sigma$ . Les variations de température entre la flamme et l'extérieur étant importantes (avec un rapport de l'ordre de 3 à 4), le coefficient de conduction $\lambda$ dépend de la température à travers une loi en puissance $\lambda(T)=\lambda_{0}T^{q}$ . Enfin pour le terme source, nous choisirons un terme de réaction constant dans la flamme et nul à l'extérieur.

En supposant que la flamme est symétrique, les conditions aux limites s'écrivent:

$$\frac{\partial T}{\partial x}(x=0)=0$$et   $$T(x=L)=T_{0}$$

$L$ étant une distance grande par rapport à l'épaisseur $\delta$ de la flamme.

En posant $u=\frac{T}{T_{0}}$ , le problème modèle associé s'écrit alors (en choisissant $L=1$ et $\delta=0,2$ ):

$$-\frac{\partial}{\partial}\left(\kappa(u)\frac{\partial u}{\partial x}\right)+\sigma(u^{4}-1)=Q(x) x\in]0,1[$$ (3.67)

$$\frac{\partial u}{\partial x}(0)=0 \„ u(1)=1$$

avec $Q(x)=\beta Heaviside(0.2-x)$ et $\kappa(u)=\kappa_{0}u^{q}$ .

3.7.2 Étude analytique

Le problème (3.67) étant non linéaire, il n'existe pas de solution analytique simple. Cependant dans le cas de faibles variations de $u$ ,

$$u\approx1 \Rightarrow\kappa(u)\approx\kappa_{0} \„ \sigma(u^{4}-1)\approx4\sigma(u-1)$$

on peut écrire un problème linéaire équivalent (en posant $\alpha=4\sigma$ ):

$$-\frac{\partial}{\partial}\left(\kappa_{0}\frac{\partial u}{\partial x}\right)+\alpha(u-1)=Q(x)$$ (3.68)

C'est le problème de diffusion linéaire avec source étudié au paragraphe c3chalsource.

Le programme Maple (3.7.2) calcule la solution générale de ce problème (3.68) en utilisant la fonction dsolve (ligne 7) et en déterminant les constantes d'intégration à l'aide des conditions aux limites (lignes 8 à 9).

programme Maple 3.7.2: solution analytique du problème 3.68

> restart;
# Diffusion dans une flamme: modele lineaire
> T0:=1;
> -diff(kappa*diff(T(x),x),x)+alpha*(T(x)-T0)=Q(x);eq:=%:
> Q:=x->beta*Heaviside(delta-x); delta:=2/10; 
# Resolution et prise en compte des C.L.
> dsolve(eq,T(x)); T1:=unapply(rhs(%),x);
> D(T1)(0)=0; T1(1)=T0; solve({%,%%},{_C1,_C2});rel:=%:
> Tex:=unapply(subs(rel,T1(x)),x);
# Trace du resultat
> subs(kappa=1/100,alpha=10,beta=30,T0=1,Tex(x));plot(%,x=0..1);

On a tracé sur la figure (3.33) cette solution analytique pour $\kappa_{0}=0,01$ et différentes valeurs de $\alpha$ pour montrer l'importance du rayonnement qui a tendance à raidir le front de température.

Figure 3.33: solution linéaire de (3.68) pour différentes valeurs de $\alpha$

3.7.3 Discrétisation par différences finies

Pour résoudre numériquement le problème non linéaire (3.67), on le discrétise tout d'abord par différences finies centrées avec un maillage de pas $dx=1/N$ sous la forme:

$$-\left(\frac{\kappa_{i+\frac{1}{2}}\left(u_{i+1}-u_{i}\right)-\ka... ...ht)}{dx^{2}}\right)+\sigma\left(u_{i}^{4}-1\right)=Q_{i} \forall i=0,N-1$$ (3.69)

en notant $\kappa_{i+\frac{1}{2}}=\kappa\frac{\kappa(u_{i+1})+\kappa(u_{i})}{2}$ et $\kappa_{i-\frac{1}{2}}=\frac{\kappa(u_{i-1})+\kappa(u_{i})}{2}$ les coefficients de diffusion en $i+\frac{1}{2}$ et $i-\frac{1}{2}$ .

A ces équations on ajoute les 2 conditions aux limites:

1. la condition de Neumann en $i=0$ , traduite comme une condition miroir pour l'équation en $i=0$ : $u_{-1}=u_{1}$ ,
2. la condition de Dirichlet en $i=N$ : $u_{N}=1$ .

On obtiens ainsi $N+1$ équations non linéaires pour $N+1$ inconnues $\{u_{i}\}_{i=0,N}$ , que l'on écrit symboliquement sous la forme:

$$F_{i}(u_{0},u_{1},..u_{j}...,u_{N})=0 \forall i=0,N$$ (3.70)

Pour résoudre numériquement ces équations, on transforme le problème en un problème équivalent de point fixe (i.e. ayant la même solution que 3.70):

$$u_{i}=G_{i}(u_{0},u_{1},..u_{j}...,u_{N}) \forall i=0,N$$ (3.71)

A partir de ce système (3.71), on construit une suite de valeurs $\{u_{i}^{k}\}_{i=0,N}$ telle que:

$$u_{i}^{k+1}=G_{i}(u_{0}^{k},u_{1}^{k},..u_{j}^{k}...,u_{N}^{k}) \forall i=1,N$$ (3.72)

Si la suite $\{u_{i}^{k}\}_{i=0,N}$ converge, elle converge vers un point fixe (i.e. une solution) de (3.71) et donc vers la solution du problème non linéaire initial (3.70). La condition de convergence de la suite de point fixe (3.72) est donnée par le théorème classique du point fixe:

théorème du point fixe:

la suite itérative (3.72) converge au voisinage d'un point fixe de (3.71) si et seulement si la matrice jacobienne $J$ de $G$ : $J_{i,j}=\left(\frac{\partial G_{i}}{\partial u_{j}}\right)$ , a une norme inférieure à 1, i.e. possède dans ce voisinage des valeurs propres de modules inférieurs à 1.

3.7.4 Schéma explicite

Une première approche consiste à considérer la solution du problème initial (3.67) comme solution stationnaire du problème instationnaire associé:

$$\frac{\partial u}{\partial t}-\frac{\partial}{\partial}\left(\kappa(u)\frac{\partial u}{\partial x}\right)+\sigma(u^{4}-1)=Q(x)$$ (3.73)

On discrétise cette équation avec du schéma différences finies explicite en temps:

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}-\left(\frac{\kappa_{i+\frac{1}{2... ...n}\right)}{dx^{2}}\right)+\sigma\left(\left(u_{i}^{n}\right)^{4}-1\right)=Q_{i}$$ (3.74)

Ce schéma explicite s'écrit s'écrit sous la forme d'un problème de point fixe:

$$u_{i}^{n+1}=G_{i}(u_{0}^{n}..u_{N}^{n})$$

avec

$$G_{i} = u_{i}^{n}+dt\left(\frac{\kappa_{i+\frac{1}{2}}^{n}\left(u_{i+1}^{... ...t)-\kappa_{i-\frac{1}{2}}^{n}\left(u_{i}^{n}-u_{i-1}^{n}\right)}{dx^{2}}\right)$$ (3.75)

$$- dt \sigma\left(\left(u_{i}^{n}\right)^{4}-1\right)+dt Q_{i}$$

La matrice jacobienne $\mathcal{G}$ de $G$ est la matrice tridiagonale d'ordre $N+1$ suivante:

$$\mathcal{G}=\left[\begin{array}{cccccc} a_{0}... ... 0 & \ddots & a_{N-1} & b_{N-1}\\ 0 & 0 & 0 & & 0 & 0\end{array}\right]$$

avec

$$a_{i} = 1-dt\left(\frac{\kappa_{i+\frac{1}{2}}^{n}+\kappa_{i-\frac{1}{2}}... ...{n}-u_{i-1}^{n}-u_{i+1}^{n})}{dx^{2}}+4\sigma\left(u_{i}^{n}\right)^{3}\right),$$

$$b_{i} = dt \frac{\kappa_{i+\frac{1}{2}+}^{n}}{dx^{2}}+\frac{\frac{1}{2}(\frac{\partial\kappa_{i+1}}{\partial u_{i+1}})(u_{i+1}^{n}-u_{i}^{n})}{dx^{2}}$$

$$c_{i} = dt \frac{\kappa_{i-\frac{1}{2}}^{n}}{dx^{2}}+\frac{\frac{1}{2}(\frac{\partial\kappa_{i-1}}{\partial u_{i-1}})(u_{i-1}^{n}-u_{i}^{n})}{dx^{2}}$$

Par définition, cette matrice Jacobienne relie une petite variation $\delta u_{i}^{n}$ de la solution à l'étape $n$ à la variation correspondante $\delta u_{i}^{n+1}$ à l'étape $n+1$ :

$$\left[\delta u_{i}^{n+1}\right]=\mathcal{G}\left[\delta u_{i}^{n}\right]$$

En interprétant $\delta u_{i}^{n}$ comme une perturbation, la relation précédente montre que la matrice $\mathcal{G}$ est la matrice d'amplification du schéma différences finies explicite (3.74), i.e. une perturbation $\epsilon_{i}^{n}$ de la solution de (3.73) vérifie

$$\left[\epsilon_{i}^{n+1}\right]=\mathcal{G}\left[\epsilon_{i}^{n}\right]$$

Or pour que le schéma (3.74) soit stable, il faut que les valeurs propres de la matrice $\mathcal{G}$ soient en module inférieures à 1. Pour notre schéma, la condition de convergence de la méthode de point fixe est en faite une condition de stabilité sur le schéma explicite associé.

Pour pouvoir faire l'analyse simplement, on simplifie les coefficients de la matrice $\mathcal{G}$ en négligeant les termes en $\frac{1}{2}(\frac{\partial\kappa_{i}}{\partial u_{i}})$ devant ceux en $\kappa_{i}$ :

$$a_{i}\approx1-dt\left(\frac{\kappa_{i+\frac{1... ...{dx^{2}}, c_{i}\approx dt \frac{\kappa_{i-\frac{1}{2}}^{n}}{dx^{2}}$$

L'équation sur la perturbation $\epsilon_{i}^{n}$ s'écrit alors:

$$\frac{\epsilon_{i}^{n+1}-\epsilon_{i}^{n}}{dt... ...silon_{i-1}^{n}\right)}{dx^{2}}+4\sigma(u_{i}^{n})^{3}\epsilon_{i}^{n}=0$$

Cette équation est analogue à l'équation de la perturbation du schéma explicite pour le problème de convection linéaire étudié au paragraphe (c3chalsource) avec les paramètres $\kappa=\kappa_{i+\frac{1}{2}}=\kappa_{i-\frac{1}{2}}$ et $\alpha=4\sigma(u_{i}^{n})^{3}$ .

L'étude de stabilité à été faîte au paragraphe 3.2.4 pour ce schéma explicite linéaire, et on a trouvé une condition de stabilité donnée par la relation (3.7).

Par analogie, la condition de stabilité pour le schéma explicite non linéaire (3.75) s'écrit donc:

$$dt<\frac{2}{4\sigma\left(u_{i}^{n}\right)^{3}+\frac{4\kappa(u_{i}^{n})}{dx^{2}}}$$

Contrairement au cas linéaire, cette condition dépend de la solution $u_{i}^{n}$ . En notant $u_{max}=Max(\vert u_{i}^{n}\vert)$ , la valeur maximale de la solution à l'étape $n$ , on choisit le pas en temps $dt$ à l'étape $n$ tel que:

$$dt<\frac{2}{4\sigma u_{max}^{3}+\frac{4\kappa(u_{max})}{dx^{2}}}$$ (3.76)

Le schéma explicite converge pour de très petites valeurs du pas en temps, vérifiant la condition (3.76). Cette condition devient très sévère lorsque le nombre de points de calcul $N$ augmente (i.e. lorsque $dx$ diminue).

3.7.5 Schéma implicite linéarisé

Pour améliorer la convergence, on peut utiliser un schéma implicite linéarisé, basé sur le schéma implicite linéaire (3.5) avec $\theta=1$ .

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt}-\left(\frac{\kappa_{i+\frac{1}{2... ...dx^{2}}\right)+\sigma\left(u_{i}^{n+1}\left(u_{i}^{n}\right)^{3}-1\right)=Q_{i}$$ (3.77)

Ce schéma s'écrit aussi sous la forme d'une itération de point fixe:

$$\left[u_{i}^{n+1}\right]=\mathcal{A}^{-1}\left[u_{i}^{n}+dt (Q_{i}+\sigma)\right]$$ (3.78)

où $\mathcal{A}$ est la matrice tridiagonale d'ordre $N+1$ suivante:

$$\mathcal{A}=\left[\begin{array}{cccccc} a_{0}... ... 0 & \ddots & a_{N-1} & b_{N-1}\\ 0 & 0 & 0 & & 0 & 1\end{array}\right]$$

avec

$$a_{i}=1+dt\left(\frac{\kappa_{i+\frac{1}{2}}^... ...2}}^{n}}{dx^{2}}, c_{i}=-dt\frac{\kappa_{i-\frac{1}{2}}^{n}}{dx^{2}}$$

La convergence de l'itération de point fixe (3.78) est lié à la stabilité du schéma (3.77). La matrice Jacobienne $\mathcal{G}$ de l'itération de point fixe s'écrit:

$$\mathcal{G}_{i,j}=\left[\mathcal{A}_{i,j}+\frac{\partial\mathcal{A}_{i,j}}{\partial u_{k}}u_{k}^{n}\right]^{-1}$$

Dans le cas linéaire, la matrice $\mathcal{A}$ est constante, et on a $\mathcal{G}=\mathcal{A}^{-1}$ . On a montré au paragraphe c3chalimp que le schéma implicite linéaire est inconditionnellement stable, et donc que la matrice $\mathcal{A}^{-1}$ possède des valeurs propres de module inférieur à 1. Dans la cas non linéaire, il faut faire intervenir la dérivée des coefficients de $\mathcal{A}$ par rapport à la solution $u_{i}^{n}$ . La matrice $\frac{\partial\mathcal{A}}{\partial u_{k}}u_{k}^{n}$ étant proportionnelle à $dt$ , pour des petits pas en temps on a $\mathcal{G\approx A}^{-1}$ , et l'itération de point fixe converge puisque les valeurs propres de $\mathcal{A}^{-1}$ sont en module plus petites que 1.

On peut donc en conclure que le schéma implicite linéarisé converge pour des petits pas en temps. Cependant la limite de convergence est beaucoup plus grande que pour le schéma explicite, et on a un schéma plus efficace.

3.7.6 Schéma de Newton

Pour rechercher les racines des équations non linéaires (3.70), on peut utiliser la méthode de Newton-Raphson, qui consiste à construite la suite itérative de point fixe suivante:

$$\left[u_{i}^{k+1}\right]=\left[u_{i}^{k}\right]-\left[\mathbf{J}^{k}\right]_{i,j}^{-1}\left[F_{j}(u_{0}^{k},..u_{N}^{k})\right]$$ (3.79)

$\mathbf{J}^{k}$ est la matrice jacobienne des fonctions $F_{i}$ : $\mathbf{J}_{i,j}^{k}=\frac{\partial F_{i}}{\partial u_{j}}(u_{0}^{k}..u_{N}^{k})$ calculée à l'itération $k$ . Cette relation s'écrit sous la forme matricielle:

$$\mathbf{J}^{k}\left[u_{i}^{k+1}-u_{i}^{k}\right]=-\left[F_{j}(u_{0}^{k},..u_{N}^{k})\right]$$ (3.80)

A chaque itération de Newton, il faut résoudre un système linéaire d'ordre $N+1$ .

Dans notre cas, la matrice Jacobienne $\mathbf{J}^{k}$ est tridiagonale et s'écrit:

$$\mathbf{J}^{k}=\left[\begin{array}{cccccc} a_... ... 0 & \ddots & a_{N-1} & b_{N-1}\\ 0 & 0 & 0 & & 0 & 1\end{array}\right]$$

avec:

$$a_{i} = \left(\frac{\kappa_{i+\frac{1}{2}}+\kappa_{i-\frac{1}{2}}}{dx^{2}... ...{n}-u_{i-1}^{n}-u_{i+1}^{n})}{dx^{2}}+4\sigma\left(u_{i}^{k}\right)^{3}\right),$$

$$b_{i} = -\frac{\kappa_{i+\frac{1}{2}}}{dx^{2}}-\frac{\frac{1}{2}(\frac{\partial\kappa_{i+1}}{\partial u_{i+1}})(u_{i+1}^{n}-u_{i}^{n})}{dx^{2}},$$

$$c_{i} = -\frac{\kappa_{i-\frac{1}{2}}}{dx^{2}}-\frac{\frac{1}{2}(\frac{\partial\kappa_{i-1}}{\partial u_{i-1}})(u_{i-1}^{n}-u_{i}^{n})}{dx^{2}}.$$

Comme pour le schéma explicite, on peut simplifier la calcul de ces coefficients en en négligeant les termes en $\frac{1}{2}(\frac{\partial\kappa_{i}}{\partial u_{i}})$ devant ceux en $\kappa_{i}$ , ce qui donne:

$$a_{i}\approx\left(\frac{\kappa_{i+\frac{1}{2}}^{n}+\kappa_{i-\fra... ...1}{2}}^{n}}{dx^{2}}, c_{i}\approx-\frac{\kappa_{i-\frac{1}{2}}^{n}}{dx^{2}}$$ (3.81)

La méthode de Newton converge au voisinage de la racine $u^{*}$ . En effet la fonction de point fixe $G_{i}$ associée s'écrit:

$$G_{i}(u_{0},..u_{N})=u_{i}-\left[\mathbf{J}\right]_{i,j}^{-1}\left[F_{j}(u_{0},..u_{N})\right]$$

et sa matrice Jacobienne $\mathcal{G}$ vaut:

$$\left[\mathcal{G}\right] = \left[I_{d}\right]-\left[\mathbf{J}\right]_{i,j}^{-1}\underbrace{... ...\left(\left[\mathbf{J}\right]\right)_{ij}^{-1}\left[F_{j}(u_{0},..u_{N})\right]$$

$$=- \frac{\partial}{\partial u_{i}}\left(\left[\mathbf{J}\right]\right)_{ij}^{-1}\left[F_{j}(u_{0},..u_{N})\right]$$

Cette matrice Jacobienne est identiquement nulle en $u=u^{*}$ car $F_{j}(u_{0}^{*},..,u_{N}^{*})=0$ . Donc sa norme vaut zéro à la racine, et la méthode de Newton converge au voisinage de cette racine.

Cette méthode de point fixe converge au voisinage de la racine (à condition que la racine soit une racine simple), et possède une vitesse de convergence quadratique. C'est la méthode la plus efficace.

3.7.7 Expérimentation numérique avec Matlab

Les trois méthodes précédentes ont été implémentées sous Matlab, et comparées sur deux cas tests:

1. cas 1: $\kappa_{0}=0.01$ $\sigma=0.1$ $\beta=1$ $\kappa(u)=\kappa_{0}\sqrt{u}$
2. cas 2: $\kappa_{0}=0.01$ $\sigma=1$ $\beta=300$ $\kappa(u)=\kappa_{0}u^{2}$

La solution calculée avec $N=51$ points est tracée sur la figure (3.34). On constate que le second cas est plus sévère que le premier, avec un profil de température plus raide.

Figure 3.34: solution du problème non-linéaire

[cas 1]

[cas 2]

3.7.7.1 Schéma explicite

Le programme Matlab (3.7.7) calcule la solution du problème non linéaire (3.67) à l'aide du schéma explicite (3.74). Le schéma est écrit pour les noeuds internes à la ligne 24, et les conditions aux limites sont imposées pour les noeuds $i=1$ (ligne 26) et $i=N$ (ligne 27). Pour tester la convergence vers la solution du problème non linéaire discret (3.70), on calcul à chaque itération en temps $n$ la norme euclidienne du résidu divisée par le nombre de points $N$ (ligne 36):

$$Err=\sqrt{\frac{1}{N}\sum_{i=1}^{N}F_{i}(u_{1}^{n},u_{2}^{n}...u_{N}^{n})^{2}}\approx\int_{0}^{1}F(u)^{2}dx$$

Cette erreur correspond à une approximation de la norme intégrale du carré du résidu $F(u)$ , qui tend vers zéro lorsque la suite converge. La solution exacte est obtenue lorsque le résidu est nul, i.e. : $F(u)=0$ . Numériquement, si ce résidu est inférieure à une valeur $\epsilon$ fixé ( $\epsilon=10^{-8}$ ), on arrête les itérations (ligne 32) et on considère que la solution numérique du problème non linéaire est la solution calculée à cette itération.

programme Matlab 3.7.7: Schèma explicite 3.74

% resolution probleme non lineaire de diffusion
% schema explicite
clear;
N=51; L=1; gamma=0.9;
dx=L/(N-1); X=[0:dx:L];
% parametre
K0=0.01; sigma=1; beta=300;
% nds internes
I=2:N-1;
% second membre
delta=0.2; Q=(X<delta)*beta;
% loi K
K=inline('t.^2');
% C.I.
Un=ones(1,N); dUn=zeros(1,N);
% iterations
nitmax=5000; epsilon=1.0e-05;
for it=1:nitmax
    % calcul du dt
    Umax=max(Un);
    dt=gamma*2/(4*sigma*Umax^3+4*K0*K(Umax)/dx^2);
    % nds internes
    K12=0.5*(K(Un(1:N-1))+K(Un(2:N)));
    dUn1(I)=K0/dx^2*(K12(I-1).*(Un(I-1)-Un(I))+K12(I).*(Un(I+1)-Un(I)))-sigma*(Un(I).^4-1.0)+Q(I);
    % C.L en 0
    dUn1(1)=K0/dx^2*(2*K12(1)*(Un(2)-Un(1)))-sigma*(Un(1).^4-1.0)+Q(1);
    dUn1(N)=0;
    % erreur
    Err=norm(dUn1)/sqrt(N);
    % fin
    Un=Un+dt*dUn1; 
    if (Err<epsilon) break; end;
end;

Le pas en temps est choisie proportionnelle à la limite de stabilité (3.76) (ligne 21).

$$dt=\gamma\frac{2}{4\sigma u_{max}^{3}+\frac{4\kappa(u_{max})}{dx^{2}}}$$ (3.82)

Pour le cas 2, nous avons testé la convergence du schéma en fonction de $\gamma$ et on a tracé le résultat sur la figure (3.35). On constate que pour $\gamma>1$ , le schéma ne converge pas $(\gamma=2)$ ou diverge $(\gamma=5)$ . Pour $\gamma<1$ , la convergence est d'autant plus grande que $\gamma$ est proche de 1. On a choisi une valeur de $\gamma\approx0.9$ proche de 1 pour avoir la meilleur convergence (ligne 4). Cela confirme l'analyse de stabilité du paragraphe 3.7.4.

Figure 3.35: convergence du schéma explicite en fonction de $\gamma$

3.7.7.2 Schéma implicite linéarisé

Le programme Matlab (3.7.7) calcule la solution du problème non linéaire (3.67) à l'aide du schéma implicite linéarisé (3.77). La structure du programme est la même que précédemment. A chaque itération en temps, on calcule le second membre et la matrice tri-diagonale du système linéaire (3.78) (lignes 29 à 32), que l'on résout à l'aide de la fonction tridiag.

programme Matlab 3.7.7: Schèma implicite linéarisé 3.77

% resolution probleme non lineaire de diffusion
% schema implicite linearise
clear;
N=51; L=1; gamma=10;
dx=L/(N-1);X=[0:dx:L];
% parametre
K0=0.01; sigma=1; beta=300;
% nds internes
I=2:N-1;
% second membre
delta=0.2; Q=(X<delta)*beta;
% loi K
K=inline('t.^2'); 
% C.I.
Un=ones(1,N); dUn=zeros(1,N); Fn=zeros(1,N);
% iterations
nitmax=1000; epsilon=1.0e-5;
for it=1:nitmax
    % calcul du dt
    Umax=max(Un);
    dt=gamma*2/(4*sigma*Umax^3+4*K0*K(Umax)/dx^2);
    % calcul du 2nd membre
    Fn(I)=Un(I)+dt*(Q(I)+sigma);
    % C.L.
    Fn(1)=Un(1)+dt*(Q(1)+sigma);
    Fn(N)=1;
    % calcul de la matrice 3Diag 
    Kn12=0.5*(K(Un(1:N-1))+K(Un(2:N))); 
    A=[1+dt*(K0/dx^2*(2*Kn12(1))+sigma*Un(1)^3),1+dt*(K0/dx^2*(Kn12(I-1)+Kn12(I))+sigma*Un(I).^3),1];
    B=[-dt*K0/dx^2*(2*Kn12(1)), -dt*K0/dx^2*(Kn12(I)),0];
    C=[0,-dt*K0/dx^2*(Kn12(I-1)),0]; 
    J=[C;A;B];
    % resolution
    Un1=tridiag(J,Fn);
    % erreur
    Err=norm((Un1-Un)/dt)/sqrt(N);
    % fin
    Un=Un1; 
    if (Err<epsilon)  break; end;
  
end;

Le pas en temps est choisie proportionnelle à la limite de stabilité (3.76) (ligne 21) suivant la relation ( ). Pour le cas 2, nous avons testé la convergence du schéma en fonction de $\gamma$ et on a tracé le résultat sur la figure (3.35).

Figure 3.36: convergence du schéma implicite en fonction de $\gamma$

On constate que contrairement au schéma implicite linéaire, le schéma n'est pas inconditionnellement stable, i.e. on ne peut pas choisir un pas en temps $dt$ arbitrairement grand, i.e. pour $\gamma>10$ , le schéma ne converge plus sans toutefois diverger. Ainsi pour $\gamma=100$ , le résidu oscille. Cette oscillation du résidu est due au choix d'un pas en temps fonction de la valeur maximum de la solution ( ). Lorsque le schéma se met à diverger, la valeur $u_{max}$ augmente et le pas en temps $dt$ diminue. Le pas en temps ayant diminué, le résidu diminue. A l'itération suivante le pas en temps augmente donc, ce qui fait croître à nouveau le résidu, et ainsi de suite. D'où la série d'oscillations sur le résidu que l'on observe sur la figure (3.35) pour $\gamma=100$ . On constate aussi que pour $\gamma\le10$ , la vitesse de convergence diminue lorsque $\gamma$ décroit. Pour $\gamma=10$ on note aussi que la vitesse de convergence est beaucoup plus grande que dans le cas du schéma explicite.

Le schéma implicite linéarisé a donc une limite de stabilité $\gamma\le10$ , mais aussi une vitesse de convergence beaucoup plus grande que le schéma explicite.

3.7.7.3 Méthode de Newton-Raphson

Le programme Matlab (3.7.7) calcule la solution du problème non linéaire (3.67) à l'aide d'une méthode de Newton (3.79). La structure du programme est la même que précédemment. Pour le calcul de la matrice Jacobienne, on a utilisé l'expression simplifiée (3.81). On note que l'utilisation d'une matrice jacobienne approchée influe uniquement sur la vitesse de convergence et non sur la solution du problème non linéaire, puisque lorsque la suite de Newton $u_{i}^{k}$ (3.80) converge elle vérifie l'équation non linéaire $F_{i}(u_{0},..u_{N})=0$ .

programme Matlab 3.7.7: Méthode de Newton-Ralphson3.77

% resolution probleme non lineaire de diffusion
% schema Newton Ralphson
clear;
N=51; L=1;
dx=L/(N-1); X=[0:dx:L];
% parametre
K0=0.01; sigma=1; beta=300;
% nds internes
I=2:N-1;
% second membre
delta=0.2; Q=(X<delta)*beta;
% loi K
K=inline('t.^2'); 
% C.I.
Un=ones(1,N); dUn=zeros(1,N); Fn=zeros(1,N);
% iterations
nitmax=1000; epsilon=1.0e-5;
for it=1:nitmax
    % calcul du 2nd membre -Fi
    % nds internes
    Kn12=0.5*(K(Un(1:N-1))+K(Un(2:N))); 
    Fn(I)=K0/dx^2*(Kn12(I-1).*(Un(I-1)-Un(I))+Kn12(I).*(Un(I+1)-Un(I)))-sigma*(Un(I).^4-1.0)+Q(I);
    % C.L en 0
    Fn(1)=K0/dx^2*(Kn12(1)*(Un(2)-Un(1))+Kn12(1)*(Un(2)-Un(1)))-sigma*(Un(1).^4-1.0)+Q(1);
    Fn(N)=0;
    % calcul de la matrice Jacobienne
    A=[K0/dx^2*(2*Kn12(1))+4*sigma*Un(1)^3,K0/dx^2*(Kn12(I-1)+Kn12(I))+4*sigma*Un(I).^3,1];
    B=[-K0/dx^2*(2*Kn12(1)), -K0/dx^2*(Kn12(I)),0];
    C=[0,-K0/dx^2*(Kn12(I-1)),0]; 
    J=[C;A;B];
    % resolution
    dUn=tridiag(J,Fn);
    % erreur
    Err=norm(Fn)/sqrt(N);
    % fin
    Un=Un+dUn;  
    if (Err<epsilon) break; end;
end;

On a comparé la vitesse de convergence de la méthode de Newton avec celle du schéma implicite linéarisé et du schéma explicite pour les deux cas tests (figure 3.37).

Figure 3.37: Comparaison de la vitesse de convergence

[cas 1]

[cas 2]

Sur cette figure, on note que la méthode de Newton est la méthode qui converge le plus rapidement, et le schéma explicite celle qui converge le plus lentement. Ainsi pour le cas 2 (le plus difficile), il faut uniquement $24$ itérations de Newton pour atteindre un résidu de $10^{-8}$ , alors qu'il en faut $316$ pour le schéma implicite linéarisé, soit $13$ fois plus, et $3552$ pour le schéma explicite, soit $148$ fois plus.