Bonjour.

Je dois résoudre une équation aux dérivées partielles (tout est dans le titre!) :
dC/dt = D * d²C/dx² - u * dC/dx
avec C=concentration
t =temps
x=distance
D= 0.01e-4 m²/s
et u = 9.910 m/s
Condition initiale : C(x, 0)= 0 à t=0.
Conditions aux limites : C (0, t)=Co = 3.5e-3 g/L et C (infini, t)=0 quelque soit t.

Pour résoudre cette équation différentielle, je souhaiterais utiliser la méthode de Crank Nicholson (schéma implicite).
Donc j’obtiens : A*Concentration (temps = t+1) = B*Concentration (temps = t).
Avec A matrice tridiagonale remplie avec les coefficients α, β et γ
B matrice tridiagonale avec les coefficients α2, β2 et γ2.

Etant donné qu’il est plus facile de faire une multiplication et une addition qu’un produit matriciel, je peux écrire B*Concentration (temps=t) comme un seul vecteur qui est la somme des termes de la concentration au temps t aux distances x+1, x, x-1.
J’ai fait attention aussi aux conditions aux limites. (C (0, t)=Co et C (infini, t)=0).

Je dois donc faire une boucle sur le temps et la distance pour trouver la concentration à t= 1 ans (= 365*24*60*60 s) et à une distance x= 36m (par exemple…) ???

Mais je n’arrive pas à résoudre ce problème avec matlab. Je ne sais pas trop où placer les boucles for. Je pense que l’une des boucles doit être comprise dans l’autre...
De plus, je ne sais pas trop comment choisir les pas de temps et de distance.

J'ai essayé de faire un programme...il me renvoie des résultats. Mais les résultats obtenus ne sont pas cohérents avec les résultats obtenus analytiquement.
Voici mon programme:

Code : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
 
N=2*36;
deltax=1; %1 mètre
deltat=24*3600; %1 jour
C0=3.5e-3;
 
%matrice A
alpha=-(1/2)*(0.01e-4/(deltax*deltax)-9.091e-7/(deltax*2))*deltat
beta=1+(0.01e-4/(deltax*deltax)*deltat)
gamma=-(1/2)*(0.01e-4/(deltax*deltax)+9.091e-7/(deltax*2))*deltat
A=zeros(N);
for i=1: (N-1)
    A(i,i)=beta;
    A(i,i+1)=alpha;
    A(i+1,i)=gamma;
end
A(N,N)=beta;
A
 
%matrice B
alpha2=(1/2)*(0.01e-4/(deltax*deltax)-9.091e-7/(deltax*2))*deltat
beta2=1-(0.01e-4/(deltax*deltax)*deltat)
gamma2=(1/2)*(0.01e-4/(deltax*deltax)+9.091e-7/(deltax*2))*deltat
B=zeros(N);
for i=1: (N-1)
    B(i,i)=beta2;
    B(i,i+1)=alpha2;
    B(i+1,i)=gamma2;
end
B(N,N)=beta2;
B;
 
%vecteur des conditions aux limites
conditionlimite=zeros(N,1);
conditionlimite(2)=-gamma*C0;
conditionlimite(N-1)=0;
conditionlimite
 
%à t=0
vectconc=zeros(N,1)
r=zeros(N,1)
for j=1:3
    vectconc=A\(vectconc+conditionlimite);
    vectconc(1)=C0;
    for i=2: (N-1)
    C=vectconc;
    r(i)=alpha2*C(i+1)+beta2*C(i)+gamma2*C(i-1);
    r(1)=C0;
    end
    vectconc=r;
end
vectconc
J’espère avoir été clair dans ma demande d’explication.
Merci beaucoup !