3.2 Équation de la chaleur avec source

3.2.1 Problème physique

On veut calculer la température $T(x,t)$ dans un barreau de longueur $L$ , de section $S$ , de masse volumique $\rho$ dont l'extrémité gauche est maintenue à température constante (i.e $T(0)=T_{e}$ ), et l'extrémité droite reçoit un flux de chaleur $\phi_{e}$ . En plus de la conduction dans le solide (le flux de chaleur par conduction dans une section s'écrit $\phi_{1}=-\lambda S\frac{\partial T}{\partial x}$ ), le barreau échange de la chaleur par convection avec l'air ambiant à la température $T_{a}$ sur toute sa longueur. En notant $\hbar$ le coefficient d'échange par convection par unité de surface, le flux de chaleur par convection s'écrit pour un élément de longueur $dx$ : $\phi_{2}=\hbar dx p (T-T_{a})$ (où $p$ est le périmètre de la section du barreau).

Figure 3.1: Température dans un barreau

L'équation de conservation de l'énergie s'écrit:

$$\rho C_{p}S\frac{\partial T}{\partial t}=\fra... ...x}\left(\lambda S\frac{\partial T}{\partial x}\right)-\hbar p (T-T_{a})$$

En notant $\kappa=\frac{\lambda}{\rho C_{p}}$ , la conductivité du matériau, $\alpha=\frac{\hbar p}{\rho C_{p}S}$ , $\Phi_{1}=\frac{\phi_{e}}{\lambda S}$ et $u=T-T_{a}$ la température relative, le problème modèle s'écrit:

Trouver $u(x,t)$ solution du système (3.1):

$$\frac{\partial u}{\partial t}=\frac{\partial}{\partial x}(\kappa\frac{\partial u}{\partial x})-\alpha u$$ (3.1)

$$u(0)=u_{0}, \frac{\partial u}{\partial x}(L)=\Phi_{1}, u(x,0)=C(x)$$

La fonction $C(x)$ correspond à la condition initiale du problème.

3.2.2 Étude de la solution analytique

Nous allons utiliser le programme Maple 3.2.2 pour déterminer la solution générale $u(x,t)$ de l'équation (3.1). Pour cela nous décomposons la solution $u(x,t)$ (lignes 4 et 5) en la somme de la solution stationnaire $u_{s}(x)$ et de la solution transitoire homogène $u_{g}(x,t)$ .

La solution stationnaire $u_{s}(x)$ est solution du problème suivant (lignes 8 à 9) (en notant $\beta=\alpha/K$ ):

$$\frac{\partial^{2}u}{\partial x^{2}}-\beta u=0$$ (3.2)

$$u(0)=u_{0}, \frac{\partial u}{\partial x}(L)=\Phi_{1}$$

que l'on résout directement avec la fonction Maple dsolve (ligne 10). Le résultat est simplifié en introduisant les fonctions hyperboliques $\sinh$ et $\cosh$ (ligne 11 à 15):

$$u_{s}(x)=\frac{\cosh(\sqrt{\frac{\alpha}{\kap... ...nh(\sqrt{\frac{\alpha}{\kappa}}x)}{\cosh(\sqrt{\frac{\alpha}{\kappa}}L)}$$

On vérifie (ligne 21) que cette solution tend vers la solution linéaire de diffusion pure lorsque $\alpha\rightarrow0$ :

$$u_{s}\longrightarrow u_{0}+\Phi_{1}x$$

programme Maple 3.2.2: Recherche de la solution analytique de l'équation 3.1

# Solution analytique de l'équation de la chaleur avec source
> restart:with(plots):
# Equation et conditions aux limites
> diff(u(t,x),t)-K*diff(u(t,x),x$2)+alpha*u(t,x)=0;eq:=%:
> u(t,0)=u0; diff(u(t,x),x)=Phi[1];
# Equation stationnaire
> beta=sqrt(alpha/K);
> -diff(U(x),x$2)+beta^2*U(x)=0; eqs:=%:
> U(0)=u0,D(U)(L)=Phi[1];CL:=%:
> dsolve({eqs,CL},U(x));rel1:=%:
> simplify(subs(exp(beta*L)=cosh(beta*L)+sinh(beta*L),
>      exp(-beta*L)=cosh(beta*L)-sinh(beta*L),
>      exp(beta*x)=cosh(beta*x)+sinh(beta*x),
>      exp(-beta*x)=cosh(beta*x)-sinh(beta*x),
> rel1));collect(%,[beta,u0]);
# solution stationnaire
> subs(beta=sqrt(alpha/K),%);sols:=%: Us:=rhs(%):
> simplify(subs(U(x)=Us,subs(beta=sqrt(alpha/K),eqs)));
> subs(K=1,alpha=1,L=1,u0=1,Phi[1]=0,Us);
> plot(%,x=0..1);
# Cas particulier alpha=0
> limit(sols,alpha=0);
# Solution instationnaire equation homogene
> diff(u(t,x),t)-K*diff(u(t,x),x$2)+alpha*u(t,x)=0;eqi:=%:
> u(t,0)=0; diff(u(t,x),x)=0;
# Separation de variables
> u(t,x)=A(t)*B(x);
> expand(subs(u(t,x)=A(t)*B(x),eqi/(A(t)*B(x))));rel2:=%:
> op(lhs(rel2))[1]+op(lhs(rel2))[3]=-op(lhs(rel2))[2]; %/K;eq1:=%:
# Les 2 membres sont constants et <0  = -omega^2
> rhs(eq1)=-omega^2; eq2:=%:
> B(0)=0,D(B)(L)=0;CL:=%:
> dsolve({eq2,B(0)=0},B(x));
> omega=(2*k+1)*Pi/2/L;rel5:=%:
> lhs(eq1)=-omega^2;%*K;
> dsolve(%,A(t));
# Solution elementaire
> sin(omega*x)*exp(-(alpha+omega^2*K)*t);
> subs(rel5,%);rel6:=%:
# Solution générale du problème transitoire
> sum(C[k]*rel6,k=0..N); Ug:=%:
# Solution exacte du probleme initiale
> Ue:=Us+Ug;

L'allure de cette solution est donnée sur la figure (3.2) (ligne 19).

Figure 3.2: Solution stationnaire du problème (3.2)

La solution transitoire $u_{g}(x,t)$ est solution du problème homogène suivant (lignes 23 et 24):

$$\frac{\partial u}{\partial t}=\kappa\frac{\partial^{2}u}{\partial x^{2}}-\alpha u$$ (3.3)

$$u(0)=0, \frac{\partial u}{\partial x}(L)=0, u(x,0)=C(x)-u_{s}(x)$$

dont on cherche la solution par une méthode de séparation de variables. On pose $u(x,t)=A(t)*B(x)$ (ligne 26) et l'équation 3.3 devient (lignes 27 à 29)

$$\frac{1}{\kappa}\left(\frac{\frac{d}{dt}A}{A}+\alpha\right)=\frac{\frac{d^{2}}{dx^{2}}B}{B}=cste$$   avec$$cste=-\omega^{2}$$

Chacun des membres de l'égalité doit être constant et négatif (on montre que si la constante est positive, alors $A=B=0$ ).

La solution $B(x)$ est obtenue directement avec la fonction dsolve (lignes 31 à 33) en utilisant la condition aux limites en $x=0$ . L'imposition de la condition en $x=L$ fixe les valeurs possibles de $\omega$ (ligne 34):

$$B(x)=B_{k}\sin\omega x$$  avec  $$\omega=\frac{2k+1}{L}\frac{\pi}{2} k=0..\infty$$

La solution $A(t)$ est aussi calculée avec dsolve (lignes 35 à 36):

$$A(t)=A_{k}e^{-(\alpha+\kappa\omega^{2})t}$$

La solution générale est une combinaison linéaire de ces solutions élémentaires (ligne 41):

$$u_{g}(x,t)=\sum_{k=0}^{\infty}C_{k}\sin\left(... ...}{L}\right) e^{-(\alpha+\kappa\left(\frac{(2k+1)\pi}{2L}\right)^{2})t}$$

d'où la solution générale $u(x,t)$ du problème (3.1)

$$u(x,t) = \frac{\cosh(\sqrt{\frac{\alpha}{\kappa}}L)\cosh(\sqrt{\frac{\alph... ...rac{\sinh(\sqrt{\frac{\alpha}{\kappa}}x)}{\cosh(\sqrt{\frac{\alpha}{\kappa}}L)}$$

$$+ \sum_{k=0}^{\infty}C_{k}\sin\left(\frac{(2k+1)\pi}{2}\frac{x}{L}\right) e^{-(\alpha+\kappa\left(\frac{(2k+1)\pi}{2L}\right)^{2})t}$$

Les coefficients $C_{k}$ sont déterminées de façon à vérifier la condition initiale $u(x,0)=C(x)$ . Ce sont les coefficients de Fourier de $C(x)-u_{s}(x)$ .

Dans la pratique, nous calculerons une solution tronquée sur $N$ modes, i.e. avec une somme finie $\sum_{k=0}^{N}$ au lieu de la somme à priori infinie $\sum_{k=0}^{\infty}$ dans l'expression précédente. Déterminons ces coefficients dans les deux cas particuliers suivants.

3.2.2.1 cas 1: solution pour $C(x)=\sin\frac{\pi x}{2L}+u_{s}$

Pour une condition initiale correspondant à la perturbation $\sin\frac{\pi x}{L}$ par rapport à la solution stationnaire, la solution pour $\alpha=1$ , $\kappa=1$ , $L=1$ , $u_{0}=1$ , $\Phi_{1}=0$ est calculée à la ligne 45 et tracé à la ligne 46 .

$$u(x,t)=\frac{\cosh(1)\cosh(x)-\sinh(1)\sinh(x)}{\cosh(1)}+\sin\frac{\pi x}{2} e^{-(1+\frac{\pi^{2}}{4})t}$$ (3.4)

En notant $t_{0}=\frac{1}{1+\pi^{2}/4}$ le temps caractéristique de diffusion, on a tracé sur la figure (3.3a) l'évolution temporelle de cette solution

programme Maple 3.2.2: Détermination de la solution analytique de l'équation 3.1 dans les cas particuliers 1 et 2

# cas particulier avec 1 seul mode
> U0:=unapply(eval(subs(N=0,Phi[1]=0,u0=1,Ue)),t,x);
> Ui:=U0(0,x);
# cas d'une condition initiale constante U=1 
> U1:=unapply(eval(subs(Phi[1]=0,u0=1,L=1,sqrt(alpha/K)=beta,Ue)),t,x);
> U1(0,x)=1; Ui:=subs(Phi[1]=0,u0=1,L=1,sqrt(alpha/K)=beta,Us);
# calcul des coefficients de Fourier Ck
> Int(Ui*sin((2*k+1)*Pi/2*x),x=0..1)+C[k]*Int(sin((2*k+1)*Pi/2*x)^2,x=0.
> .1)=Int(sin((2*k+1)*Pi/2*x),x=0..1);
# Calcul des differents termes de cette expression
> k:='k':int(Ui*sin((2*k+1)*Pi/2*x),x=0..1);
> simplify(subs(sin(k*Pi)=0,%));rel8:=%:
> simplify(numer(rel8)/expand(denom(rel8)));rel9:=%:
> 2*(2*k+1)*Pi/(4*beta^2+Pi^2*(2*k+1)^2);rel10:=%:
> simplify(rel9-rel10);
> int(sin((2*k+1)*Pi/2*x),x=0..1);simplify(subs(sin(k*Pi)=0,%));
> rel11:=%:
> int(sin((2*k+1)*Pi/2*x)^2,x=0..1);simplify(subs(sin(k*Pi)=0,%));
> rel12:=%:
> rel12*C[k]=rel11-rel10;
# D'ou l'expression des Ck et d la solution 
> C[k]=(rel11-rel10)/rel12; rel13:=rhs(%):
> U2:=unapply(subs(C[k]=rel13,beta=1,K=1,N=20,U1(t,x)),t,x);

Figure 3.3: Évolution temporelle de $u(x,t)$ pour $\kappa =1$ , $\alpha =1$ , $\Phi _{1}=0$

[cas 1][cas 2]

3.2.2.2 cas 2: solution pour $C(x)=1$

Pour déterminer la solution, on calcule les coefficients $C_{k}$ en multipliant la condition initiale $u_{s}(x)+u_{g}(x,0)=C(x)$ par $\sin\left(\frac{2k+1}{2L}\pi x\right)$ et en intégrant de $x=0$ à $L$ (lignes 51 à 53). On laisse Maple faire les intégrations (lignes 55 à 64) en simplifiant les expressions et on obtiens les coefficients $C_{k}$ (lignes 67):

$$C_{k}=\frac{4}{(2k+1)\pi}-\frac{4(2k+1)\pi}{4\frac{\alpha^{2}}{\kappa^{2}}+\pi^{2}(2k+1)^{2}}$$

On trace l'allure de la solution calculée avec $N=20$ modes sur la figure (3.3b) pour $t$ variant de 0 à $2t_{0}$ .

3.2.3 Discrétisation avec un $\theta$ -schéma

Le $\theta$ -schéma est obtenue par combinaison linéaire d'un schéma explicite ( $\theta=0$ ) et d'un schéma implicite ( $\theta=1$ ) . Il s'écrit pour l'équation (3.1):

$$\frac{u_{i}^{n+1}-u_{i}^{n}}{dt} = (1-\theta)\left(\kappa\frac{u_{i+1}^{n}-2u_{i}^{n}+u_{i-1}^{n}}{dx^{2}}-\alpha u_{i}^{n}\right)$$

$$\theta\left(\kappa\frac{u_{i+1}^{n+1}-2u_{i}^{n+1}+u_{i-1}^{n+1}}{dx^{2}}-\alpha u_{i}^{n+1}\right)$$ (3.5)

Le $\theta$ -schéma est un schéma implicite (sauf pour $\theta=0$ ) qui couple les valeurs inconnues $\{u_{i}^{n+1},u_{i+1}^{n+1},u_{i-1}^{n+1}\}$ à l'étape $t^{n+1}$ aux valeurs connues $\{u_{i}^{n},u_{i+1}^{n},u_{i-1}^{n}\}$ à l'étape $t^{n}$ comme indiqué sur le diagramme (3.4).

Figure 3.4: diagramme du $\theta$ -schéma pour l'équation de la chaleur

Nous allons tout d'abord étudier les propriétés de ce schéma en utilisant Maple.

3.2.4 Étude de la stabilité du $\theta$ -schéma

Dans le programme Maple (3.2.4), on effectue l'étude de la stabilité du $\theta$ -schéma (3.5). Pour cela on définit l'équation exacte (ligne 3), et l'équation discrétisée (ligne 6), que l'on stocke dans les variables eq et eqh. La perturbation $U_{p}(n,i)$ est définie à la ligne 10, comme le mode de Fourier $\psi_{n}e^{I\omega idx}$ . L'équation étant homogène, l'équation sur la perturbation est l'équation discrète eqh dans laquelle on remplace la solution numérique $U[n,i]$ par la perturbation $U_{p}(n,i)$ (ligne 11). On simplifie cette équation (ligne 12), et l'on en déduit le facteur d'amplification $G=\frac{\psi_{n+1}}{\Psi_{n}}$ . Pour simplifier l'expression de $G$ , on pose $y=\cos(\omega dx)$ et on introduit les 2 paramètres $r$ et $s$ suivants:

$$r=\kappa\frac{dt}{dx^{2}} \„ s=\alpha dt$$ (3.6)

On obtiens ainsi le facteur d'amplification $G$ (ligne 17):

$$G=y\rightarrow\frac{2 yr\theta-2 yr-1+2 r-2 r\theta+s-s\theta}{2 yr\theta-1-2 r\theta-s\theta}$$

La condition de stabilité impose:

$$\left\vert G(y)\right\vert\le1 \forall y\in[-1,+1]$$

programme Maple 3.2.4: Etude de la stabilité du $\theta$ -schéma 3.5

# Etude du theta-schema pour l'équation de la chaleur
> restart:
> eq:=diff(U(t,x),t)-kappa*diff(U(t,x),x$2)+alpha*U(t,x)=0:eq;
# Discretisation DF (theta schema)
> delta2:=(f,n)->(f[n,i+1]-2*f[n,i]+f[n,i-1])/dx^2;
> eqh:=(U[n+1,i]-U[n,i])/dt-kappa*(theta*delta2(U,n+1)+(1-theta)*
   delta2(U,n))+theta*alpha*U[n+1,i]+(1-theta)*alpha*U[n,i]=0:
> eqh;
# Etude de la stabilite: introduction d'une perturbation
> Up:=(n,i)->Psi[n]*exp(I*omega*i*dx);
> subs(U[n,i]=Up(n,i),U[n,i+1]=Up(n,i+1),U[n,i-1]=Up(n,i-1),
    U[n+1,i]=Up(n+1,i),U[n+1,i+1]=Up(n+1,i+1),
    U[n+1,i-1]=Up(n+1,i-1),eqh);
> eqhp:=%:
> expand(simplify((eqhp)*dt*exp(-I*omega*i*dx)));rel1:=%:
> G:=solve(simplify(subs(Psi[n+1]=G*Psi[n],rel1/Psi[n])),G);
> G:=simplify(subs(cos(omega*dx)=y,alpha=s/dt,kappa=r*dx^2/dt,G));
# Schema explicite  theta=0
> G0:=unapply(subs(theta=0,G),y);
> plot(subs(r=1/2,s=1/2,G0(y)),y=-1..1);
> G0(1)>-1;
> G0(-1)>-1;
> s<2;2>s+4*r;cdt1:=%:
> subs(r=kappa*dt/dx^2,s=alpha*dt,cdt1);
# Condition de stabilité pour le schema explicite
> (alpha+4*kappa/dx^2)*dt<2;
# Schema implicite theta=1
> G1:=unapply(subs(theta=1,G),y);
> plot(subs(r=1,s=1,G1(y)),y=-1..1);
> simplify(G1(-1));simplify(G1(1));
# Schema inconditionellement stable
# Schema C.N. theta=0.5
> G2:=unapply(subs(theta=1/2,G),y);
> plot(subs(r=1,s=1,G2(y)),y=-1..1);
> simplify(G2(-1));simplify(G2(1));
> G2(1)>-1;numer(G2(1))>-denom(G2(1));
> G2(-1)>-1; numer(G2(-1))>-denom(G2(-1));
# Schema inconditionellement stable

Pour simplifier, nous ferrons l'analyse pour les 3 valeurs suivantes de $\theta$ : $\theta=0, \frac{1}{2}, 1$ .

3.2.4.1 $\theta=0$ schéma explicite

Dans ce cas, l'expression de $G$ est donnée à la ligne 19

$$G=y\rightarrow2 yr+1-2r-s$$

C'est une fonction croissante de $y$ , inférieure à 1 pour $y\in[-1,+1]$ . La condition $\left\vert G\right\vert\le1$ implique donc $G(-1)>-1$ , qui s'écrit (ligne 23):

$$s+4r\leq2$$

En exprimant $s$ et $r$ en fonction de $dt$ et $dx$ , on obtiens la condition de stabilité (ligne 24):

$$\left(\alpha+\frac{4\kappa}{dx^{2}}\right)dt\leq2$$

Cette condition impose une valeur maximum du pas en temps $dt$ pour un maillage en espace donnée (i.e. pour $dx=\frac{L}{N}$ fixé):

$$dt\leq\frac{2}{\alpha+\frac{4\kappa}{dx^{2}}}$$ (3.7)

Pour $\alpha=0$ , on retrouve la condition de stabilité classique $\kappa\frac{dt}{dx^{2}}\le\frac{1}{2}$ . Pour $\kappa=0$ , on a une équation différentielle en temps et on retrouve la condition de stabilité de la méthode d'intégration d'Euler $\alpha dt\leq2$ .

Le schéma explicite est donc conditionnellement stable avec une condition de stabilité donnée par la relation (3.7).

3.2.4.2 $\theta=1$ schéma implicite

Pour $\theta=1$ , l'expression de $G$ est donnée à la ligne 28

$$G=y\rightarrow\frac{1}{1+2 r (1-y)+s}$$

C'est une fonction croissante de $y$ , qui est comprise entre 0 et 1 pour $y\in[-1,+1]$ . On le vérifie facilement en traçant la fonction $G$ (ligne 29) et en calculant $G(-1)$ et $G(1)$ (ligne 30) .

Le schéma implicite est donc inconditionnellement stable.

3.2.4.3 $\theta=\frac{1}{2}$ schéma de Crank Nicholson

Pour $\theta=\frac{1}{2}$ , l'expression de $G$ est donnée à la ligne 33

$$G=y\rightarrow\frac{1-r(1-y)-\frac{s}{2}}{1+r(1-y)+\frac{s}{2}}$$

C'est une fonction croissante de $y$ , qui est comprise entre -1 et 1 pour $y\in[-1,+1]$ . On le vérifie facilement en traçant la fonction $G$ (ligne 34) et en calculant $G(-1)$ et $G(1)$ (lignes 35 à 37) .

Le schéma de Crank Nicholson est donc inconditionnellement stable.

3.2.5 Étude de la consistance du $\theta$ -schéma

Pour étudier la consistance du $\theta$ -schéma, nous utiliserons le programme Maple (3.2.5). Pour cela on définit l'équation exacte (ligne 3), et l'équation discrétisée (ligne 6), que l'on stocke dans les variables eq et eqh.

programme Maple 3.2.5: Etude de la consistance du $\theta$ -schéma 3.5

# Etude du theta-schema pour l'équation de la chaleur
> restart:
> eq:=diff(U(t,x),t)-kappa*diff(U(t,x),x$2)+alpha*U(t,x)=0:eq;
# Discretisation DF (theta schema)
> delta2:=(f,n)->(f[n,i+1]-2*f[n,i]+f[n,i-1])/dx^2;
> eqh:=(U[n+1,i]-U[n,i])/dt-kappa*(theta*delta2(U,n+1)+(1-theta)*
   delta2(U,n))+theta*alpha*U[n+1,i]+(1-theta)*alpha*U[n,i]=0:
> eqh;
# Etude de la consistence
# calcul de l'erreur de troncature
> Uex:=(p,q)->U(t+(p-n)*dt,x+(q-i)*dx);
> subs(U[n,i]=Uex(n,i),U[n,i+1]=Uex(n,i+1),U[n,i-1]=Uex(n,i-1),
   U[n+1,i]=Uex(n+1,i),U[n+1,i+1]=Uex(n+1,i+1),
   U[n+1,i-1]=Uex(n+1,i-1),lhs(eqh));
> rel3:=%:
> expand(simplify(rel3-lhs(eq)));rel4:=%:
# Developpement en serie de Taylor
> U(t+dt,x)=convert(mtaylor(U(t+dt,x),[dt],4),diff);S1:=%:
> U(t+dt,x+dx)=convert(mtaylor(U(t+dt,x+dx),[dt,dx],4),diff);
  S2:=%:
> U(t+dt,x-dx)=convert(mtaylor(U(t+dt,x-dx),[dt,dx],4),diff);
  S3:=%:
> U(t,x+dx)=convert(mtaylor(U(t,x+dx),[dx],5),diff);S4:=%:
> U(t,x-dx)=convert(mtaylor(U(t,x-dx),[dx],5),diff);S5:=%:
> simplify(subs(S1,S2,S3,S4,S5,rel4)):
> ErrT:=collect(%,dt);
# Calcul du coefficient du terme en dt
> coeff(ErrT,dt);c1:=%:
> diff(eq,t);
> factor(simplify(c1+theta*lhs(%)));
# Donc schema en O(dt,dx^2) si theta#1/2 et O(dt^2,dx^2) sinon

On remplace ensuite dans l'équation discrétisée eqh la solution numérique $U[n,i]$ par la solution exacte $U_{ex}(n,i)$ définie à la ligne 11. La fonction $U_{ex}$ est définie de sorte à faire apparaître les solutions à $t+dt$ , $x-dx$ ou $x+dx$ . Cette substitution est faite à la ligne 12 et le résultat est stocké dans la variable rel3. On calcul l'erreur de troncature en soustrayant à rel3 l'équation exacte eq et on stocke le résultat dans rel4. On effectue ensuite des développements en série de Taylor autour de la solution $u(x,t)$ (lignes 18 à 24), que l'on substitue dans rel4 (ligne 25). L'erreur de troncature est obtenue en l'écrivant sous forme d'un polynôme en $dt$ (ligne 26):

$$ErrT = 1/6 \alpha \theta \frac{\partial^{3}}{\partial t^{3}}U(t,x) \textit{dt}^{3}$$

$$+ \left(1/6 \frac{\partial^{3}}{\partial t^{3}}U(t,x)+1/2 \alpha \theta \frac{\partial^{2}}{\partial t^{2}}U(t,x)\right)\textit{dt}^{2}$$

$$+ \left(1/2 \frac{\partial^{2}}{\partial t^{2}}U(t,x)-\kappa \the... ...al x}U(t,x)+\alpha \theta \frac{\partial}{\partial t}U(t,x)\right)\textit{dt}$$

$$- 1/12 \kappa \left(\frac{\partial^{4}}{\partial x^{4}}U(t,x)\rig... ... \theta \left(\frac{\partial^{4}}{\partial x^{4}}U(t,x)\right)\textit{dx}^{2}$$

On remarque que le coefficient du terme en $dt$ ressemble à la dérivée temporelle de l'équation exacte. On va donc le transformer en utilisant l'équation exacte. On extrait ce coefficient de la variable ErrT (ligne 28). On dérive l'équation exacte par rapport à t (ligne 29) et on la soustraie du coefficient après multiplication par $\theta$ (ligne 30). Le coefficient du terme en $dt$ s'écrit alors (ligne 20):

$$\frac{1}{2}(2\theta-1)\frac{\partial^{2}}{\partial t^{2}}U(t,x)$$

Ce coefficient s'annule pour $\theta=\frac{1}{2}$ .

En conclusion le $\theta$ -schéma (3.5) est consistant avec l'équation (3.1) et l'erreur de troncature est en $O(dt,dx^{2})$ sauf pour $\theta=\frac{1}{2}$ (schéma de Crank Nicolson) pour lequel l'erreur de troncature est en $O(dt^{2},dx^{2})$ .

3.2.6 Expérimentation numérique avec Matlab

Le $\theta$ -schéma s'écrit sous la forme suivante suivante (avec les paramètres $r$ et $s$ donnés en (3.6)):

$$-\theta ru_{i-1}^{n+1}+(1+\theta(s+2r))u_{i}^{n+1}-\theta ru_{i+1}^{n+1} = (1-\theta)ru_{i-1}^{n}$$ (3.8)

$$+ (1-(1-\theta)(s+2r))u_{i}^{n}$$

$$+ (1-\theta)ru_{i+1}^{n}$$

Cette équation n'est évidemment valable que pour les noeuds internes $i=1,N-1$ . Pour les noeuds $i=0$ et $i=N$ il faut prendre en compte les conditions aux limites. En $x=0$ ( $i=0$ ), on a une condition de Dirichlet, et donc

$$u_{0}^{n+1}=u_{0}$$ (3.9)

En $x=L$ ( $i=N$ ), on a une condition de Neumann, qu'il faut discrétiser pour calculer $u_{N}^{n+1}$ . Pour cela on utilise une condition miroir, qui consiste à calculer une approximation centrée de la dérivée au point $i=N$ , pour approximer la valeur manquante $u_{N+1}$ à l'aide de la condition aux limites:

$$\left.\frac{\partial u}{\partial x}\right\ver... ...N+1}-u_{N-1}}{2 dx} \Longrightarrow u_{N+1}=u_{N-1}+2\Phi_{1}dx$$

Cela nous permet d'écrire l'équation (3.8) pour $i=N$ sous la forme:

$$-2\theta ru_{N-1}^{n+1}+(1+\theta(s+2r))u_{N}^{n+1} = 2(1-\theta)ru_{N-1}^{n}$$ (3.10)

$$+ (1-(1-\theta)(s+2r))u_{N}^{n}$$

$$+ 2r\Phi_{1}dx$$

Les équations (3.8 à 3.10) peuvent s'écrire sous la forme d'un système matriciel:

$$\mathcal{A}\left[\begin{array}{c} u_{0}^{n+1} u_{1}^{n+1} u_{... ...}\right]+\left[\begin{array}{c} 0 0 0 \vdots 0 \Phi\end{array}\right]$$ (3.11)

où $\mathcal{A}$ et $\mathcal{C}$ sont les 2 matrices tridiagonales suivantes (en notant $\theta_{1}=1-\theta$ et $\Phi=2r\Phi_{1}dx$ ):

$$\mathcal{A} = \left[\begin{array}{cccccc} 1 & 0 & 0 & & 0 & 0\\ -\theta r & 1+... ...+2r) & -\theta r\\ 0 & 0 & 0 & & -2\theta r & 1+\theta(s+2r)\end{array}\right]$$

$$\mathcal{C} = \left[\begin{array}{cccccc} 1 & 0 & 0 & & 0 & 0\\ \theta_{1}r & ... ...theta_{1}r\\ 0 & 0 & 0 & & 2\theta_{1}r & 1-\theta_{1}(s+2r)\end{array}\right]$$

Nous avons programmé ce $\theta$ -schéma (3.11), en utilisant une programmation Matlab optimisée. C'est la fonction tshema (programme 3.2.6), qui calcule la solution numérique dans le vecteur T pour une condition initiale T $_{i}$ , des valeurs $\theta$ , $r$ , $s$ et $\Phi$ données et un nombre d'itérations en temps nit fixé. Elle utilise un algorithme spécifique aux systèmes tri-diagonaux: l'algorithme de Thomas (décrit dans l'annexe chap3D). C'est un algorithme classique pour la résolution de systèmes tri-diagonaux, mais il n'est pas implémenté dans Matlab. Nous avons donc écrit une fonction en language C tridiag avec une interface Matlab qui est décrite au paragraphe lin3D.

programme Matlab 3.2.6: fonction tshema: utilisation de l'algorithme de Thomas

function [T]=tshema(Ti,theta,r,s,phi,nit)
% calcul de la solution T du theta-schema 
% Ti=solution initiale
% theta,r,s,phi parametres du shema
% nit= nbre de pas d'integration
% methode de Thomas avec tridiag
N=length(Ti);
% matrices 3 Diag
e=ones(1,N);
A=[-theta*r*e; (1+theta*(2*r+s))*e; -theta*r*e];
C=[(1-theta)*r*e; (1-(1-theta)*(2*r+s))*e; (1-theta)*r*e];
% C.L dirichlet en 1, Neuman en N (cdt miroir)
A(:,1)=0; A(2,1)=1;          C(:,1)=0; C(2,1)=1;
A(1,N)=2*A(1,N); A(3,N)=0;   C(1,N)=2*C(1,N); C(3,N)=0;
% solution initiale
T=Ti;
% iterations
for it=1:nit
 B=C(1,1:N)'.*[0;T(1:N-1)]+C(2,1:N)'.*T(1:N)+C(3,1:N)'.*[T(2:N);0];
 B(N)=B(N)+phi;
 T=tridiag(A,B);
end

La fonction Matlab tshema (3.2.6) utilise cette fonction tridiag. Pour cela on construit les matrices tridiagonales $\mathcal{A}$ et $\mathcal{C}$ comme des matrices rectangulaires $3*N$ (lignes 10 à 11). On notera aussi la forme du produit matriciel $\mathcal{C}.T^{n}$ (ligne 19) pour le calcul du second membre $B$ .

En utilisant cette structure de données et l'algorithme de Thomas tridiag, on obtiens une méthode très efficace, qui pour $N=1000$ est $250$ fois plus rapide que la méthode de Gauss classique implémentée dans Matlab. On trouvera dans l'annexe chap3D, une comparaison de l'efficacité des différentes méthodes de résolution de systèmes tridiagonaux disponibles sous Matlab par rapport à cet algorithme de Thomas.

Avec les conditions initiales des cas 1 et 2, $\kappa=1$ , $\alpha=1$ , $\Phi_{1}=0$ , on a calculé la solution numérique avec $\theta=1/2$ , un maillage de $N=50$ points et $r=0.1$ . Le résultat est tracé sur la figure (3.5). En comparant avec les solutions exactes (figures 3.3), on constate que l'on a calculé une bonne approximation de la solution.

Figure 3.5: Solution numérique avec $N=50$ , $r=0.2$ , $\theta =\frac {1}{2}$ pour $\kappa =1$ , $\alpha =1$

[cas 1]

[cas 2]

Pour quantifier l'erreur entre la solution exacte et la solution numérique, nous allons faire une étude numérique de la précision du $\theta$ -schéma pour l'équation de la chaleur (3.1).

3.2.7 Étude numérique de la précision

Pour étudier numériquement la précision du $\theta$ -schéma (3.11) , nous avons effectuer deux séries de calcul: une pour déterminer la précision spatiale et une pour déterminer la précision temporelle.

3.2.7.1 Précision spatiale du $\theta$ -schéma

On a comparé la solution numérique avec la solution analytique (3.4) au bout d'un temps $t=t_{0}/2$ ( $t_{0}$ étant le temps caractéristique de diffusion du problème $\frac{4}{\pi^{2}}$ ). On a choisi un pas d'intégration $\Delta t$ suffisamment petit pour que l'erreur d'intégration en temps soit négligeable par rapport à l'erreur spatiale (on a fait le calcul avec une valeur de $r$ fixée et égale à $0.1$ ). On a tracé sur la figure (3.6a) l'erreur maximale en fonction de $N$ pour les 3 valeurs $\theta=0,\frac{1}{2},1$ . On vérifie sur la figure que l'erreur est d'ordre 2 pour tous les schémas (i.e. en $O(N^{-2})$ ).

Figure 3.6: précision spatiale et temporelle du $\theta$ -schéma (3.11)

[précision spatiale]

[précision temporelle]

3.2.7.2 Précision temporelle du $\theta$ -schéma

L'étude de la précision temporelle est un peu plus délicate, car il faut choisir un nombre de points $N$ et des pas en temps $\Delta t$ tels que l'erreur globale du schéma soit contrôlée par l'erreur d'intégration en temps. Cela nécessite des calculs avec des pas en temps grands, qui ne vérifient plus la condition de stabilité du schéma explicite. Nous avons donc fait l'étude uniquement pour $\theta=1$ et $\theta=\frac{1}{2}$ , en comparant la solution analytique (3.4) au bout d'un temps $t=1.26 t_{0}$ avec la solution calculée numériquement. Pour les valeurs de $\Delta t$ choisies, le paramètre $r$ varie de $256$ à $1$ , et est donc bien en deçà de la limite de stabilité du schéma explicite $(\theta=0)$ . On a tracé sur la figure (3.6b) l'erreur maximale en fonction de $\Delta t$ . On vérifie sur la figure que le schéma implicite $(\theta=1)$ est d'ordre 1 en temps (i.e. en $O(\Delta t)$ ), et que pour les grands $\Delta t$ le schéma Cranck Nicholson $(\theta=\frac{1}{2})$ est d'ordre 2 en temps (i.e. en $O(\Delta t^{2})$ ). Pour les petits $\Delta t$ , on note une saturation de l'erreur du schéma Cranck Nicholson, qui devient alors contrôlée uniquement par la précision spatiale.