Résolution d'équation de convection pure

Salut à tout ! j'aimerais résoudre en 1D & 2D l'équation :
dC/dt+VdC/dx=0 C :inconnu & V :vitesse de transport
en servant le shéma centré ,shéma amont,shéma LAX-Wendroff,shéma LAX-Wendroff avec limiteur.
Merci !

heu.... sans vouloir te vexé, ta question ne relève pas de fortran....

commence par discrétiser ton équation via les schémas que tu as cité sur papier (c'est des maths pas de la programmation) ensuite tu passe au codage et franchement, l'avantage du fortran c'est que t'as (persque) qu'à recopié tes equations....

Salut

Justement ça c'est mon projet en simulation par F77 et que j'ai eu beaucoup de problèmes car le sujet ne te permet pas de comprendre avec exactitude de ton problème !
Par exemple pour le shéma centré on a
C(i,n+1)=C(i,n)-Sigma*V*[C(i+1,n)-C(i-1,n)]
avec Sigma=Deltat/2Deltax
Deltax=XL/(Nx-1)
les données sont : XL :taille du domaine
Nt :nb de pts dépent du temps t
Nx :nb de pts dépent du position x
Deltat=0.01
V(i)=V0 :vitesse de trasport est constant
Condition initiale : f(x)=0.5+0.5*tanh(80*(x-0.1)) <=> C(i,0)=f(x(i))
Pour la CI j'ai fais :
Code:

1 2 3 4 5 do i=1,Nx V(i)=V0 x(i)=[(2*i-1)*Deltax]/2 C(i)=0.5+0.5*tanh(80*(x(i)-00.1)) end do
-- cette CI ne dépent pas de temps pour tant C dépent de x et de t
-- par la suite j'aimerais bien compresser en une seul dimension en faisant
Code:

1 2 3 4 5 do p=1,Nt do q=1,Nx C(q)=C(q)-Sigma*V0*[C(q+1)-C(q-1)] end do end do
--là j'ai un problème de borne , donc il faut que je fasse
Code:

1 2 3 4 5 6 7 do p=1,Nt C(1)=C(0)-Sigma*V0*[C(2)-C(0)] -> ??? C(Nx)=C(q)-Sigma*V0*[C(q+1)-C(q-1)] ->??? do q=2,Nx-1 C(q)-Sigma*V0*[C(q+1)-C(q-1)] end do end do
--donc ma question est : quelle technique permet-on compresser en 1D ? et comment résoudre le problème de borne ? rapellons que on veut avoir :

Citation:

Nt

1 C1
.
.
i Ci
.
.
Nx C(Nx)

bon pour commencer, arrete le fortran 77 c'est un peu dinosaure!!! (tu veras le format libre de fortran ça change la vie!!!! :D ) cela dit:

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
 
program convection
implicit none
integer,parameter :: ki=kind(1.d0)   !dimension des real comme des double precision
integer,parameter :: NT=10000 !dimension des tableaux en temp
integer,parameter :: NX=250 !dimension en x
 
real(ki),parameter :: XL=10.d0        !longueur de l'échantillon
real(ki),parameter :: DeltaT=1.d-2   !pas de temps
real(ki),parameter :: DeltaX=XL/(Nx+1) !pas d'espace
real(ki),paramater :: sigma=Deltat/(2.d0*Deltax)
real(ki),paramater :: V0=10.d0        !vitesse initiale (que tu peux rentrer après mais pour simplifié je la pose au départ....
real(ki),dimension(0:NX,0:NT) :: C        !notre inconnu fonction de l'espace et du temps
 
integer i,nn  ! compteurs en espace et en temps
 
!=====================================================
! début du code
!=====================================================
!initialisation
 
init:do i=0,NX
  C(i,0)=5.d-1*(1.d0+dtanh(8.d1*(i*DeltaX-1.d-1)))
end do init
 
!condition de bords: (ici je prend une condition de diriclet homogène car tu n'a pas renseigné ce que sont les bords....)
 
bord:do nn=0,NT
  C(0,nn)=0.d0
  C(Nx,nn)=0.d0
end do bord
 
! boucle temporelle (je te fait un schémas d'euler explicit mais tu peut 
!l'améliauré rapidement en runge-kutta classique d'ordre 4 avec pas de temps 
!adaptatif si tu veux)
Boucletemp:do nn=0,NT-1
  bclSpace:do i=1,Nx-1
    C(i,nn+1)=C(i,nn)+sigma*V0*(C(i+1,nn)-C(i-1,nn))
  end do bclSpace
end do Boucletemp
 
!====================================================
! stockage
! là je te fait un stockage en version matricielle avec les x en horizontal et le t en vertical:
 
open(unit=20,file='C.dat',status='unknown')
write(20,*)'C '                                          !on dit ce qi'il y a dans le fichier
write(20,1000)' x=',(i*dx,i=0,Nx)                      !on écrit l'échelle en x
write(20,*)'t='                                          !décorations
bcltempo:do nn=0,NT                                !ligne par lignes
  write(20,1001)(nn*DeltaT),(C(i,nn),i=0,Nx)     ! les valeur de C
end do bcltempo
close(20)                                                !jamais oublier de fermé le fichier
1000  format(a3,256d25.16)                       !les format utilisés....
1001  format(2500e25.16)
 
END program convection

pour résoudre ton truc il te faut des condition des bords, on ne résout pas un problème aux limites sans connaitre les limites....

Merci Genteur Slayer de m'avoir donné des indications pourtant ce n'est pas vraiment ce que tu as compris !
--Pour le constant C,il est plutot conseille de faire un vecteur 1D C(0:Nx) et d'ecraser à chaque pas de temps les valeurs de C(i) ! on sauvegarde les valeurs de C a quelques instants, par exemple toutes les 100 iterations temporelles
--Pour le condition de bords on utilise un schema decentre (on utilise que les points à l'interieur du domaine) et on l'integre a la boucle de temps, ecrive le developpement limite de C(x+dx) et C(x-dx) et tu vas voir.Et
Code:

1 2 3 4 do nn=0,Nt C(0,nn)=0 ! aucune raison que ce soit 0 C(Nx,nn)=0 ! aucune raison que ce soit 0 end do
-- Boucle temporelle
Code:

1 2 3 4 5 6 7 8 9 10 11 12 do nn=0,Nt-1 ! ici on fait en i=0, en utilisant les valeurs en i=0 et i=1 ! C(0)= C(0) + 2*sigma*V*(C(1)-C(0)) ! on multiplie sigma par 2 car on est entre i=0 et i+1=1, donc on a Deltax et non pas 2*Deltax au denominateur do i=1,Nx-1 C(i,nn+1)=C(i,nn)+sigma*V*(C(i+1,nn)-C(i-1,nn)) end do ! ici on fait en i=Nx, en utilisant les valeurs en i=Nx-1 et i=Nx ! C(Nx)= C(Nx) + 2*sigma*V*(C(Nx)-C(Nx-1)) ! stockage des valeurs, si nn multiplie de 100 end do
--Stockage :Je dois représenter graphiquement les données dans un ficher par itération avec en colonne x, (y en 2D), C(x,y)
8O

faire comme tu dis pour C peut parfois être dangereux (dans ce cas précis non) car tu a parfois plusieurs variables interdépendantes et donc tu as besoin à tout moment des prédécésseurs temporels (exemple, méca flotte) cependant dans ton cas c'est pareil, juste tu économise de la mémoire (environ nx*(nt-1)*4 o soit 9.536Mo dans le cas que je t'es donné (à condition que tu compile sur un pc standard)

de plus pour utiliser le shéma avant en début de x et le shémas arrirère en fin de x, tu n'a qu'as le faire en dehors de la boucle central
Code:

1 2 3 4 5 6 7 8 9 10 ! boucle temporelle (je te fait un schémas d'euler explicit mais tu peut !l'améliauré rapidement en runge-kutta classique d'ordre 4 avec pas de temps !adaptatif si tu veux) Boucletemp:do nn=0,NT-1 C(0,nn+1)=C(0,nn)+2.d0*sigma*V0*(C(1,nn)-C(0,nn)) !<<< schem av en debut x C(Nx,nn+1)=C(Nx,nn)+2.d0*sigma*V0*(C(Nx,nn)-C(Nx-1,nn)) !<<< schem arr en fin de x bclSpace:do i=1,Nx-1 C(i,nn+1)=C(i,nn)+sigma*V0*(C(i+1,nn)-C(i-1,nn)) end do bclSpace end do Boucletemp
en virrant la partie juste avant cette boucle dans ce que je t'es donné pour gagner quelques milisecondes de calcul.

je vois pas où est le soucis

à moins que tu ne cherche pas un schémas explicit en temps

En fait, ce que je voudrais fais c'est de bien savoir comment est ce possible une fonction qui dépend du temps & d'espace peut être compressé en 1D et voilà la solution que j'ai trouvé :

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
 
      program projet
      implicit none
      external f1
      real f1
      integer nx,nt,i
      real tf,deltat,deltax,pi,x,t,c,ti,pf,sigma,v
      dimension x(0:10000),t(0:10000),c(0:10000)
      t(0)=0
      ti=t(0)
 
      write(*,*) 'choisir une vitesse initiale V'
      read(*,*) v
      write(*,*) 'choisir le temps final'
      read(*,*) tf
      write(*,*) 'choisir la position initiale'
      read(*,*) pi
      write(*,*) 'choisir la position finale'
      read(*,*) pf
      write(*,*) 'choisir le nombre de pts de position'
      read(*,*) nx
      deltax=(pf-pi)/nx
      write(*,*) 'choisir le nombre de pts de temps'
      read(*,*) nt
 
      deltat=(tf-ti)/nt
      sigma=deltax/(2*deltat)
      x(0)=pi
      do i=0,nx
      x(i+1)=x(i)+deltax
      end do
      do i=0,nt
      t(i+1)=t(i)+deltat
      end do
      call centre(f1,nx,nt,deltax,v,c,sigma)
      end
 
      subroutine centre(f,nx,nt,deltax,v,c,sigma)
      dimension c(0:10000)
      integer nt,nx,i
      real c,deltax,f,sigma,v
 
      do i=0,nx
      c(i)=f(i*deltax)
      end do
 
      do j=1,nt
      c(0)=c(0)-2*sigma*v*(c(1)-c(0))
      do i=1,nx-1
      c(i)=c(i)-sigma*v*(c(i+1)-c(i-1))
      end do
      c(nx)=c(nx)-2*sigma*v*(c(nx)-c(nx-1))
      end do
 
      do i=0,nx
      write(*,*) 'c(',i,')= ',c(i)
      end do
      return
      end
 
      real function f1(x)
      real x
      f1=0.5+0.5*tanh(80*(x-0.1))
      return
      end

! Le temps ici c'est le temps pour d'effectuer les calculs donc avec ce code ça marche ! Que pense -tu de ça ? Ah oui autrefois j'ai trompé de singe de c(0)=c(0)+2*sigma*v*(c(1)-c(0)) au lieu de c(0)=c(0)-2*sigma*v*(c(1)-c(0)) mais ça ne change rien puis qu'on n'intéresse que la partie de code , non ? :aie:

J'ai oublié de faire une graphe maic car je n'ai pas fini la moitié de mon projet donc je vais faire plustard !:mouarf:

je vois un petit problème:
Code:

1 2 3 4 5 6 7 8 do j=1,nt c'(0)=c(0)-2*sigma*v*(c(1)-c(0)) do i=1,nx-1 c'(i)=c(i)-sigma*v*(c(i+1)-c(i-1)) end do c'(nx)=c(nx)-2*sigma*v*(c(nx)-c(nx-1)) end do
c'est ton code, mais j'ai mis des primes (') sur les variable à t+1 j'arrive pas à dire ça plus clairement: quand tu rentre dans ta boucle do (en i): à la première itération tu as :
c'(1)=c(1)-sigma*v*(c(1)-c(0)) or si on regarde la ligne du dessus, c(0)=c'(0) du coup tu as des valeur qui ne provienne pas le même itération temporelle: ce que je te conseil pour resté dans cette méthode c'est d'avoir deux tableaux c et de les échanger une fois le calcul terminé:
Code:

1 2 3 4 5 6 7 8 9 10 do j=1,nt c2(0)=c(0)-2*sigma*v*(c(1)-c(0)) do i=1,nx-1 c2(i)=c(i)-sigma*v*(c(i+1)-c(i-1)) end do c2(nx)=c(nx)-2*sigma*v*(c(nx)-c(nx-1)) !là je le rajoute c=c2 !je sais pas si ça marche en F77 car je crois que c du F95 (copie globale de tableau) end do
au pire tu fait une boucle:
Code:

1 2 3 4 do alpha=0,nx c(alpha)=c2(alpha) end do
de plus tu n'as pas du tout besoin des tableau x et t, il ne te serve (si j'ai bien compris que pour le stockage) tu peux les faire calculer sur place (c'est pas très lourd et comme ça tu économise la ram) même si cette optimisation est une affaire de gout... à toi de voir

ensuite je comprend pas ce que tu veux dire par 'compressé en 1D' tu est déjà en 1D (instationnaire certes)....

pour ton stockage, c'est toi qui choisit ton format: tu peut écrire dans un fichier un truc style:
Code:

1 2 3 4 5 6 7 8 9 t= 1 2 3 4 5 ......... x= 1 c(1,1) c(1,2) c(1,3 ..... 2 c(2,1) .... 3 ... 4 5 ...
ou bien un truc genre (g écris c en deux variable mais c juste pour la compréhention)
Code:

1 2 3 4 5 6 7 8 9 10 t= x= c(x,t)= 1 1 c(1,1) 1 2 c(2,1) 1 ... ... ... ... ... 2 1 c(1,2) 2 2 c(2,2) 2 ... .... 2
c'est toi qui choisit le format en fonction de ce que tu as besoin en sortie cependant avec ton code si tu veux tous les (ou quelques) pas de temps il te faudra faire des stockages DANS la boucle en j si tu choisit de faire c en deux variables, tu peux faire tous tes calculs et après stocker
si tu n'a besoin que de la dernière itération en temps alors: choisit le second format de fichier (tu peux même le stocké sur une ligne au lieu de la colonne) et stocke après le calcul....

je refais encore mon prg et je vais te fournir l'énoncé

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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
program projet
      implicit none
      external f1
      integer nx,nt,i,choix,choix1
      real*8 tf,deltat,deltax,pi,x,t,c,ti,pf,sigma,v,f1,gamma
      dimension x(0:10000),t(0:10000),c(0:10000),gamma(0:10000)
      parameter (pi=0)
12    write(*,*)'choissir une vitesse initiale v'
      read(*,*,err=12) v
18    write(*,*)'choissir la position finale'
      read(*,*,err=18) pf
20    write(*,*)'choissir le nombre de pts de position'
      read(*,*,err=20) nx
22    write(*,*)'choissir le nombre de pas de temps'
      read(*,*,err=22) nt
24    write(*,*)'choissir une condition initiale : ',
     &          'CI 1 taper 1',
     &          'CI 2 taper 2'
      read(*,*,err=24) choix1
26    write(*,*)'choissir un shema :                                ',
     &          'shema centre : taper 1                             ', 
     &          'shema amont : taper 2                              ',
     &          'shema Lax_Wendroff:  taper 3', 
     &          'shema Lax_Wendroff avec le limiteur de flux : taper 4'
      read(*,*,err=26) choix     
 
      deltax=(pf-pi)/nx        !XL=f-pi
      deltat=deltax/2*v        !boite noie
      sigma=deltax/(2*deltat) !
 
      open(10,file='projet.dat',status='old')
 
      do i=0,nt
      t(i)=i*deltat
      write(10,*) 't(',i,')=',t(i)
      end do
 
      do i=0,nx
      x(i)=i*deltax
      write(10,*) 'x(',i,')= ',x(i)
      end do
 
      if (choix.eq.1) then
                call centre(f1,nx,nt,deltax,v,c,sigma)
      elseif (choix.eq.2) then
                call amont(f1,nx,nt,deltax,v,c,sigma)
      elseif (choix.eq.3) then
                call Lax_Wendroff1(f1,nx,nt,deltax,v,c,sigma)
      elseif(choix.eq.4) then
               call Lax_Wendroff2(f1,nx,nt,deltax,v,c,sigma,gamma)
      else
         go to 26
      end if
      end
 
      subroutine centre(f1,nx,nt,deltax,v,c,sigma)
      external f1
      dimension c(0:10000),c_old(0:10000)
      integer nt,nx,i,alpha
      real*8 c,c_old,deltax,f1,sigma,v
 
      do i=0,nx
      c_old(i)=f1(i*deltax)
      end do
 
      do j=0,nt
         do i=1,nx-1
         c(i)=c_old(i)-sigma*v*(c_old(i+1)-c_old(i-1))
         end do
         c(0)=c(1)
         c(nx)=c(nx-1)
         c_old(i)=c(i)
      end do
 
      do i=0,nx
         write(10,*) 'c(',i,')= ',c(i)
         write(*,*) 'c(',i,')= ',c(i)
      end do
      return
      end
 
      subroutine amont(f1,nx,nt,deltax,v,c,sigma)
      external f1
      dimension c(0:10000),c_old(0:10000)
      integer nx,nt,i,j
      real*8 q,c,deltax,sigma,v,f1,c_old
      parameter (q=0)                                
 
       do i=0,nx
       c_old(i)=f1(i*deltax)
       end do
 
       do j=0,nt
       do i=1,nx-1
       c(i)=c_old(i)-sigma*(max(q,v)*(c_old(i)-c_old(i-1))
     &  + min(q,v)*(c_old(i+1)-c_old(i)))
       end do
       c(0)=c(1)
       c(nx)=c(nx-1)
       c_old(i)=(i)
       end do
 
       do i=0,nx
       write(10,*) 'c(',i,')= ',c(i)
       write(*,*) 'c(',i,')= ',c(i)
       end do
       return
       end
 
      subroutine Lax_Wendroff1(f1,nx,nt,deltax,v,c,sigma)
      external f1
      dimension c(0:100000),c_old(0:100000)
      integer nt,nx,i,j
      real*8 c,c_old,deltax,f1,sigma,v
 
       do i=0,nx
          c_old(i)=f1(i*deltax)
       end do
 
       do i=0,nt
       do j=1,nx-1
       c(i)=c(i)-sigma*v*(c_old(i)-c_old(i-1))
     &   +1/2*sigma*v*(v*sigma-1)*(c_old(i+1)-2*c_old(i)-c_old(i-1))
       end do
       c(0)=c(1)
       c(nx)=c(nx-1)
       c_old(i)=c(i)
       end do
 
       do i=0,nx
          write(10,*) 'c(',i,')= ',c(i)
          write(*,*) 'c(',i,')= ',c(i)
       end do
       return
       end 
 
      subroutine Lax_Wendroff2(f1,nx,nt,deltax,v,c,sigma,gamma)
      dimension c(0:100000),c_old(0:100000),gamma(0:100000)
      external f1
      integer nt,nx,i,j
      real*8 c,c_old,deltax,f1,sigma,q,gamma,p,v
      parameter (p=0,q=1)
 
      do i=0,nx
        c_old(i)=f1(i*deltax)
      end do
 
      do j=0,nt
         gamma(0)=0
         do i=1,nx-1
         gamma(i)=max(q,min(p,2*((c_old(i)-c_old(i-1))
     &   /(c_old(i+1)-c_old(i)))))*(c_old(i+1)-c_old(i))/deltax
         end do
         c(0)=c(1)
         c(nx)=c(nx-1)
         c_old(i)=c(i)
         gamma(nx)=0
      end do
 
      do j=0,nt
         do i=1,nx-1
         c(i)=c_old(i)-sigma*max(q,v)*(c_old(i)-c_old(i-1))
     &   +1/2*sigma*deltax*max(q,v)*(v*sigma-1)*(gamma(i)-gamma(i-1))
       end do
         c(0)=c(1)
         c(nx)=c(nx-1)
         c_old(i)=c(i)
       end do
 
       do i=0,nx
       write(10,*) 'c(',i,')= ',c(i)
       write(*,*) 'c(',i,')= ',c(i)
       end do
       return
       end
 
!Condition initiale
 
      real*8 function f1(x)
      real*8 x
      if (choix1.eq.1) then
            f1=0.5+0.5*tanh(80*(x-0.1))
      elseif (choix1.eq.2) then
        if ((x .ge. 0.) .and. (x .lt. 0.2)) then
            f1 = 1
         elseif ((x .ge. 0.2) .and. (x .lt. 0.6)) then
            f1 = 0
         elseif ((x .ge. 0.6) .and. (x .lt. 0.8)) then
            f1 = sin(pi*(x-0.6)/0.2)**2
         elseif (x .ge. 0.8) then
           f1 = 0
         end if
      end if
        return
      end

Et voici l'énoncé:
1ere partie:
http://rapidshare.com/files/29618627/doc1.jpg
2ieme partie:
http://rapidshare.com/files/29618453/doc2.jpg
3ieme partie:
http://rapidshare.com/files/29618653/doc3.jpg

slt ! c'est bon j'ai réussi de le faire car je suis en période d'exam donc je ne peux pas poster ma correction alors attends moi quelques jours :yaisse2:

Comme promise voilà mon program

Citation:

program PROJET2(1D)
external f1,f21,f22,f23,f24
real XL,deltax,c,f1,v,deltat,sigma,x,cold
integer nx,nt,choix,choix1
dimension c(0:10000),x(0:10000),cold(0:10000)

1 write(*,*)'choissir un shema : ',
& 'shema centre : taper 1 ',
& 'shema amont : taper 2 ',
& 'shema Lax_Wendroff : taper 3 ',
& 'shema Lax_Wendroff avec le limiteur de flux : taper 4'
read(*,*,err=1) choix1

2 write(*,*)'choissir une condition initiale : ',
& 'CI 1 taper 1',
& 'CI 2 taper 2'
read(*,*,err=2) choix

3 write(*,*)'choissir une vitesse initiale v'
read(*,*,err=3) v ! v=v0 : vitesse constante
4 write(*,*)'choissir la longueur XL'
read(*,*,err=4) XL
5 write(*,*)'choissir le nombre de pts de position'
read(*,*,err=5) nx
6 write(*,*)'choissir le nombre de pas de temps'
read(*,*,err=6) nt

!L'enregistrement des différents shémas aux conditions initiales proposées

open(1,file='centre1.txt',status='old')
open(2,file='centre2.txt',status='old')
open(3,file='amont1.txt',status='old')
open(4,file='amont2.txt',status='old')
open(5,file='lw1.txt',status='old')
open(6,file='lw2.txt',status='old')
open(7,file='lwf1.txt',status='old')
open(8,file='lwf2.txt',status='old')

deltax=XL/nx
deltat=deltax/(2*v) !critère de stabilité
sigma=deltat/(2*deltax)
! Paramètres de la condition initiale 2
x1=0.2
x2=0.6
x3=0.8

do i=0,nx
x(i)=i*deltax
end do
! la comparaison des choix
select case(choix)
case(1)
if (choix1.eq.1) then
call centre1(f1,c,deltax,nx,nt,sigma,v,cold)
elseif (choix1.eq.2) then
call amont1(f1,c,deltax,nx,nt,sigma,v,cold)
elseif (choix1.eq.3) then
call lw1(f1,c,deltax,nx,nt,sigma,v,cold)
elseif(choix1.eq.4) then
call lwlf1(f1,c,deltax,nx,nt,sigma,v,cold)
else
go to 1
end if
case(2)
do i=1,nx
if (x(i).ge.0.and.x(i).lt.x1) then
c(i)=f21(i*deltax)
elseif (x(i).ge.x1.and.x(i).lt.x2) then
c(i)=f22(i*deltax)
elseif (x(i).ge.x2.and.x(i).lt.x3) then
c(i)=f23(i*deltax)
elseif (x(i).ge.x3) then
c(i)=f24(i*deltax)
end if
end do
if (choix1.eq.1) then
call centre2(nx,nt,v,c,sigma,cold)
elseif (choix1.eq.2) then
call amont2(nx,nt,v,c,sigma,cold)
elseif (choix1.eq.3) then
call lw2(nx,nt,v,c,sigma,deltax,cold)
elseif(choix1.eq.4) then
call lwlf2(c,deltax,nx,nt,sigma,v,cold)
else
go to 1
end if
case default
go to 2
end select
end

! Schéma centré à condition initiale 1

subroutine centre1(f1,c,deltax,nx,nt,sigma,v,cold)
external f1
real c,f1,deltax,cold
dimension c(0:10000), cold(0:10000)!cold est un tableau auxilliaire
! qui permettra de stocker les valeurs de C temporellement
! Calcul de C initial à partir de f1
do i=0,nx
c(i)=f1(i*deltax)
cold(i)=c(i)
end do

do j=1,nt !boucle temporelle dont on exclue la CI déjà calculée
do i=1,nx-1
c(i)=c(i)-sigma*v*(c(i+1)-cold(i-1))
end do

do i=0,nx !réinitialisation de cold afin de pouvoir le faire
cold(i)=c(i) !évoluer en même temps que C et pour qu'il ne se fixe à
end do !une seule valeur
c(0)=c(1) !conditions de Neuman
c(nx)=c(nx-1)
end do

do i=0,nx
write(1,*) c(i)
end do

return
end

! Schéma centré à condition initiale 2

subroutine centre2(nx,nt,v,c,sigma,cold)
dimension c(0:10000), cold(0:10000)
integer nt,nx,i,j
real cold,c,sigma,v

do i=0,nx
cold(i)=c(i)
end do

do j=1,nt
do i=1,nx-1
c(i)=c(i)-sigma*v*(c(i+1)-cold(i-1))
end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(2,*) c(i)
end do

return
end

! Schéma Amont à condition initiale 1

subroutine amont1(f1,c,deltax,nx,nt,sigma,v,cold)
external f1
dimension c(0:10000),cold(0:10000)
real c,f1,deltax,cold

do i=0,nx
c(i)=f1(i*deltax)
cold(i)=c(i)
end do

do j=1,nt
do i=1,nx-1
c(i)=c(i)-2*sigma*(max(0.,v)*(c(i)-cold(i-1))+min(0.,v)*
$(c(i+1)-c(i)))
end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(3,*) c(i)
end do

return
end

! Schéma Amont à condition initiale 2

subroutine amont2(nx,nt,v,c,sigma,cold)
dimension c(0:10000),cold(0:10000)
integer nx,nt,i,j
real c,v,sigma,cold

do i=0,nx
cold(i)=c(i)
end do

do j=1,nt
do i=1,nx-1
c(i)=c(i)-2*sigma*(max(0.,V)*(c(i)-cold(i-1))+ min(0.,V)*
$ (c(i+1)-c(i)))
end do
c(0)=c(1)
c(nx)=c(nx-1)

do i=0,nx
cold(i)=c(i)
end do

end do

do i=0,nx
write(4,*) c(i)
end do

return
end

! Schéma Lax-Wendroff à condition initiale 1

subroutine lw1(f1,c,deltax,nx,nt,sigma,v,cold)
external f1
dimension c(0:10000),gamma(0:10000),cold(0:10000)
real c,f1,deltax,gamma,b,cold

do i=0,nx
c(i)=f1(i*deltax)
cold(i)=c(i)
end do

do j=1,nt
gamma(0)=0 ! initialisation de gamma nécessaire afin de ne pas
do i=1,nx-1 !engendrer des erreurs
gamma(i)=(c(i+1)-c(i))/deltax
b=gamma(i)-gamma(i-1)
gamma(i-1)=(c(i)-cold(i-1))/deltax
c(i)=c(i)-2*sigma*v*(c(i)-cold(i-1))
$ +0.5*2*sigma*deltax*(V*sigma-1)*b
end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(5,*) c(i)
end do

return
end

! Schéma Lax-Wendroff à condition initiale 2

subroutine lw2(nx,nt,v,c,sigma,deltax,cold)
dimension c(0:10000),gamma(0:10000),cold(0:10000)
integer nx,nt,i,j
real c,v,sigma,b,gamma,cold

do i=0,nx
cold(i)=c(i)
end do

do j=1,nt

gamma(0)=0
do i=1,nx-1
gamma(i)=(c(i+1)-c(i))/deltax
gamma(i-1)=(c(i)-cold(i-1))/deltax
b=gamma(i)-gamma(i-1)
c(i)=c(i)-2*sigma*V*(c(i)-cold(i-1))+sigma*V*deltax*(V*sigma-1)*b
end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(6,*) c(i)
end do

return
end

! Schéma Lax-Wendroff avec limiteur de flux à condition initiale 1

subroutine lwlf1(f1,c,deltax,nx,nt,sigma,v,cold)
external f1
dimension c(0:10000),cold(0:10000)
real c,f1,deltax,g1,g2,tt1,tt2,tt3,g3

dx=deltax
s=2*sigma

do i=0,nx
c(i)=f1(i*deltax)
cold(i)=c(i)
end do

do j=1,nt
do i=1,nx-1

tt1=(c(i)-cold(i-1))/(c(i+1)-c(i))
tt2=(cold(i-1)-cold(i-2))/(c(i)-cold(i))
tt3=(c(i+1)-c(i))/(C(i+2)-c(i+1))

g1=max(0.,min(1.,2*tt1))*(c(i+1)-c(i))/dx
g2=max(0.,min(1.,2*tt2))*(c(i)-cold(i))/dx
g3=max(0.,min(1.,2*tt3))*(c(i+2)-c(i+1))/dx

c(i)=c(i)-s*max(0.,v)*(c(i)-cold(i))+0.5*s*dx*max(0.,v)*(v*s-1)*
$(g1-g2)-s*min(0.,v)*(c(i+1)-c(i))+0.5*s*dx*min(0.,v)*(v*s+1)*
$(g3-g1)

end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(7,*) c(i)
end do

return
end

! Schéma Lax Wendroff avec limiteur de flux à condition initiale 2

subroutine lwlf2(c,deltax,nx,nt,sigma,v,cold)
real c,deltax,g1,g2,tt1,tt2,tt3,g3,cold
dimension c(0:10000), cold(0:10000)

dx=deltax
s=2*sigma

do i=0,nx
cold(i)=c(i)
end do

do j=1,nt
do i=1,nx-1

tt1=(c(i)-cold(i-1))/(c(i+1)-c(i))
tt2=(cold(i-1)-cold(i-2))/(c(i)-cold(i))
tt3=(c(i+1)-c(i))/(c(i+2)-c(i+1))

g1=max(0.,min(1.,2*tt1))*(c(i+1)-c(i))/dx
g2=max(0.,min(1.,2*tt2))*(c(i)-cold(i))/dx
g3=max(0.,min(1.,2*tt3))*(c(i+2)-c(i+1))/dx

c(i)=c(i)-s*max(0.,v)*(c(i)-cold(i))+0.5*s*dx*max(0.,v)*(v*s-1)*
$(g1-g2)-s*min(0.,v)*(c(i+1)-c(i))+0.5*s*dx*min(0.,v)*(v*s+1)*
$(g3-g1)

end do

do i=0,nx
cold(i)=c(i)
end do
c(0)=c(1)
c(nx)=c(nx-1)
end do

do i=0,nx
write(8,*) c(i)
end do
return
end

! Les functions ultilisées

real function f1(x)
f1=0.5+0.5*tanh(80*(x-0.1))
return
end
real function f21(x)
f21=1
return
end
real function f22(x)
f22=0
return
end
real function f23(x)
f23=(sin((3.14*(x-0.6))/0.2))**2
return
end
real function f24(x)
f24=0
return
end

Pour le graphe j'ai une petite prolème car je l'ai travaillé avec EXEL 2007 et qu'ici (salle informatique d'université ) on ultilise EXEL 2003.donc désolé pour les graphes !

juste pour dire que pour tracer des graph, Excel y a bcp mieux et même dans les gratuit (gnuplot, medit...) et surtout je te conseil de sortir tous tes documents en pdf....
moi, g reçu autant de format de fichier que d'élèves dans mes classes et c'est profondément chiant de trouver quel soft ouvre ces fichiers... au moins, avec le pdf t'es sûr que y aura pas de soucis de mise en page, version ou quoi....
avec pdfcreator, qui est gratuit, tu sort tes fichier word, excel, openoffice... en pdf...

enrgeticien

Citation:

Envoyé par truong

Salut à tout ! j'aimerais résoudre en 1D & 2D l'équation :
dC/dt+VdC/dx=0 C :inconnu & V :vitesse de transport
en servant le shéma centré ,shéma amont,shéma LAX-Wendroff,shéma LAX-Wendroff avec limiteur.
Merci !

j'ai votre programme en 1D
avec shéma