[Débutant] equation de diffusion avec matlab [Résolu]

[VauRDeC]

Bonjour,
Je cherche à résoudre l'équation de diffusion (de la température) dans un champs 2D plus précisément avec les différence finies (sous matlab)...
Si vous avez autres chose je suis preneur également...
Merciiii a+

[Dut]

Tu devrais commencer par nous montrer ce que tu as déjà codé sous MATLAB (même si c'est faux) et nous dire clairement le point qui te bloque.

Nous ne sommes en aucun la pour résoudre cet exercice à ta place.

[VauRDeC]

Envoyé par Dut

Nous ne sommes en aucun la pour résoudre cet exercice à ta place.

En fait je cherchai juste si vous auriez quelque exemples de résolution qui pourraient me servir de support, c'est tout et c'est clair que je ne cherche pas à ce que vous fassiez mon travail...

Alors voila, c'est une équation de diffusion instationnaire :dTanaly/dt - ( d2Tanaly/dx2 + d²Tanaly/dy² ) = Surs(x,y,t)..
Tanaly= T analytique = T(x,y,t)=sin((x*k/a))*sin((y*kp/b))*cos(o*t)
Surs= Source

Les dérivées secondes spatiales sont approchées par les différence centrée à l'ordre, et la différentiation temporelle par un schéma d'Euler retardé..
Les conditions aux limites sont de types Dirichlet..

Voici le code:

Code :

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
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
%Programme equation de diffusion stationnaire 2D
clear all
z=1;
w=1;
%définition du maillage a N+1 points entre o et 1
N=7;
M=7;
x=0:z/(N-1):z;
dx=z/(N-1);
y=0:w/(M-1):w;
dy=w/(M-1);
dt=0.005;
a=2/(dx*dx)+2/(dy*dy)+3/(2*dt);
b=-1/(dx*dx);
c=-1/(dx*dx);
d=-1/(dy*dy);
e=-1/(dy*dy);
f=1/dt;
o=2*3.145159/f;
k=1;
kp=1;
fi=0;
fip=0;

%défintion de l'opérateur

rft=125
operat=zeros(N*M);
operat(1:M,1:M)=eye(M);
operat(N*M-M+1:N*M,N*M-M+1:N*M)=eye(M);

% Définition des sous-pérateurs opDiag, opbeta, opgama

% Définition des sous-pérateurs opDiag, opbeta, opgama

opDiag=eye(M);
for i=1
    opDiag(i,i)=0;
end
for i=M
    opDiag(i,i)=0;
end

opbeta=zeros(M);
for i=1:M-1
    opbeta(i,i+1)=1;
end
for i=1
    opbeta(i,i+1)=0
end

opgama=zeros(M);
for i=1:M-1
    opgama(i+1,i)=1;
end
for i=M-1
    opgama(i+1,i)=0;
end

% Définition des matrices pour les sous-opérateurs gama, beta, Diag
%sous matrice gama
gama1=eye(M);
for i=1
    gama1(i,i)=0;
end
for i=M
    gama1(i,i)=0;
end
gama=gama1*c;
%sous matrice beta
beta1=eye(M);
for i=1
    beta1(i,i)=0;
end
for i=M
    beta1(i,i)=0;
end
beta=b*beta1;
%pour la definition de la sous matrice Diag, on l'écrit comme la somme de
%trois matrice: alphadM, epsdM et deltadM
alphadM=a*eye(M);
for i=1
    alphadM(i,i)=1;
end
for i=M
    alphadM(i,i)=1;
end

epsdM=zeros(M);
for i=1:M-1
    epsdM(i,i+1)=e;
end
for i=1
    epsdM(i,i+1)=0;
end

deltadM=zeros(M);
for i=1:M-1
    deltadM(i+1,i)=d;
end
for i=M-1
    deltadM(i+1,i)=0;
end
deltadM=deltadM;


Diag=alphadM+epsdM+deltadM;

masterop=operat+kron(opDiag,Diag)+kron(opbeta,beta)+kron(opgama,gama);

%fin de définition de notre opérateur.

sdg=M*N

%Initialisation du champ Tnm1 au temps initial t=0 sur tout le domaine

v=0
for i=1:N;
    for j=1:M
        u=(i-1)*M+j;
        Tnm1(u)=sin(k*i+fi)*sin(kp*j+fip)*cos(o*dt*v);
    end
end

%Initialisation du champ Tn au temps t=dt sur tout le domaine

v=1
for i=1:N
    for j=1:M
	u=(i-1)*M+j;
        Tn(u)=sin(k*i+fi)*sin(kp*j+fip)*cos(o*dt*v)
    end
end

% Calcul pour le 3eme champ Tnp1 au pas de temps suivant t=2dt

for v=2:10

for i=1:N;
    for j=1:M;
        u=(i-1)*M+j;
        Tanaly(u)=sin(k*i)*sin(kp*j)*cos(o*v*dt));
        Surs(u)=sin(k*i)*sin(kp*j)*(cos(o*v*dt)0(k*k+kp*kp)-o*sin(o*v*dt));
        secondmembr(u)=S(u)+4/(2*dt)*Tn(u)-1/(2*dt)*Tnm1(u);
    end
end

% Les conditions aux limites

for j=1:M;
    i=1;
    u=(i-1)*M+j;
    secondmembr(u)=Tanaly(u);
    i=N;
    u=(i-1)*M+j;
    secondmembr(u)=Tanaly(u);
end

for i=2:N-1
    j=1;
    u=(i-1)*M+j;
    secondmembr(u)=Tanaly(u);
    j=M;
    u=(i-1)*M+j;
    secondmembr(u)=Tanaly(u);
end
 

%Résolution de l'équation
 
 Minv=inv(masterop);
  Tnp1=(inv(masterop)*secondmembr')';
Tnm1=Tn;
Tn=TNp1;

end

En fait je bloque sur l'itération du champ temporelle, ça ne marche pas ou si ça marche j'ai une très grosse erreur sur mon T cherché.

PS:La définition de l'opérateur (masterop) est très longue, je sais il y a plus simple mais j'ai choisi de définir plusieurs sous-matrice et ainsi de déterminer mon opérateur par le produit tensoriel
Si vous avez des question sur mon raisonnement ou que j'ai oublié de donné des infos je suis toute ouïe

MErciiii
Bne journée a+