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
| clear all
clf
clf
%equ du/dt-D d²u/dx²=0;
%homogène en x=0
%chauffe en Rb u(t)=tanh(t) qui varie entre 0 & 1
%condition initiale homogène
%Inputs
choix=input('Voulez un schéma explicite (0) ou implicite (défaut)?');
Lb=-10; %left boundary
Rb=10; %right boundary
Lx=Rb-Lb; %Largeur domaine
N=100; %nombre d'éléments géométriques
n=N+1; %nombre de noeuds
pas=(Rb-Lb)/N; %taille d'un élément spatial h
Nt=1000; %nombre de noeuds en temps
Lt=1; %taille de l'intervalle temporel
Dt=Lt/(Nt-1); %pas de temps
V=100; %vitesse
%maillage spatial
Vmaill=[Lb:pas:Rb];
%Condition initiale
for i=1:n
U0(i)=(1+cos(2*pi/Lx*((i-1)*pas+Lx/2)))/2; %condition initiale
end
a=V*Dt/2/pas %condition de stabilité en explicite
if choix==0
figure(1)
hold on
plot(Vmaill,U0,'m')
%Schéma numérique
%Explicite
%Bords (Périodicité)
AE(1,1)=1;
AE(1,2)=-a;
AE(1,n)=+a;
AE(n,1)=-a;
AE(n,n-1)=a;
AE(n,n)=1;
%intérieur
%matrice AE
for i=2:N
AE(i,i)=1 ;
AE(i,i-1)=+a;
AE(i,i+1)=-a;
end
%U(temps+1)=A*U(temps)
VE=U0'; %vecteur de travail explicite
for k = 1:Nt+1
VV=AE*VE; %ittération temporelle
VE=VV; %stockage résultat
if rem(k,10)==0
figure(1)
hold on
plot(Vmaill,VE,'b')
end
end
else
figure(2)
hold on
plot(Vmaill,U0,'m')
%Schéma Implicite
%Bords (Périodicité)
AI(1,1)=1;
AI(1,2)=a;
AI(1,n)=-a;
AI(n,1)=a;
AI(n,n-1)=-a;
AI(n,n)=1;
%intérieur
%matrice AI
for i=2:N
AI(i,i)=1 ;
AI(i,i-1)=-a;
AI(i,i+1)=a;
end
%U(temps+1)=A*U(temps)
VE=U0'; %vecteur de travail explicite
for k = 1:Nt+1
VV=inv(AI)*VE; %ittération temporelle
VE=VV; %stockage résultat
if rem(k,10)==0
figure(2)
hold on
plot(Vmaill,VE,'b')
end
end
end
disp('FIN DU PROGRAMME') |
Partager