résolution d'un systémé linéaire par la méthode de Thomas

**adadu26** · 28/10/2018, 23h00

Bonjour,

J'ai fais un code Matlab pour résoudre un problème d'advection de l'air grâce aux différences finies. Mon code est le suivant. J'ai cependant un problème, je ne comprend pas pourquoi ma variable X vaut toujours 0?? Est ce un problème sur ma boucle for ou autre? je n'arrive pas a trouver et je voulais savoir votre avis.

Merci de l'attention accordé a monmessage

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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
clear all 
 
 
%Vitesse de propagation
V=1;
 
%Longueur du domaine
L=2;
 
%définition du nombre de mailles
Nx=200;
 
%calcul de la taille des mailles
dX=L/(Nx-1);
 
%vecteur pour la representation graphique
vecX=zeros(Nx,1);
for i=1:Nx
    vecX(i)=(i-1)*dX;
end
 
%pas de temps
dt=0.01;
 
% Initialisation du problème
C_old=zeros(Nx,1);
for i=1:Nx
    C_old(i)=exp(-5*(5*(i-1)*dX-1)^2);
end 
%plot(vecX,C_old);
hold on
 
%Boucle sur les pas de temps
tic
for t=1:80
    %remplissage de la matrice z et du vecteur B avec les bonnes expressions
    z=zeros(Nx,3);
    B = zeros(Nx,1);
    for i=1:Nx
        if (i==1) %condition en amont
            z(i,2)=1;
            B(i)=0;
        elseif (i==Nx) %condition en aval
            z(i,2)=1;
            B(i)=0;
        else  %Equation au coeur du domaine
            z(i,2)=1/(dt);
            z(i,3)=V/(4*dX);
            z(i,1)=-V/(4*dX);
            B(i)=C_old(i)/(dt)-V/(4*dX)*(C_old(i+1)-C_old(i-1));
        end
    end
 
%Resolution du Système Linéaire
%[A,B]=F_GaussPiv(A,B);
%[C_old]=Tr_Sup(A,B);
 
alpha=zeros(Nx,1);
beta=zeros(Nx,1);
 
alpha(1)=z(1,2);
for i=2:Nx
    beta(i)=z(i,1)/alpha(i-1);
    alpha(i)=z(i,2)-beta(i)*z(i-1,3);
end
 
 
y=zeros(Nx,1);
y(1)=B(1);
for i=2:Nx
    y(i)=B(i)-beta(i)*y(i-1);
end
 
X=zeros(Nx,1);
X(Nx)=y(Nx)/alpha(Nx);
 
for i=Nx-1:1
    X(i)=(y(i)-z(i,3)*X(i+1))/alpha(i);
end
 
 
 
 
%Stockage des valeurs au pas de temps précédent
if mod(t,20)==0
    plot(vecX,X)
    hold on
end
end
 
 
toc