Discussion :

[bengir]

Bonjour,
J'arrive tout juste sur le forum donc je m'excuse d'avance si je ne procède pas comme il faut.
Voilà mon problème, je dois résoudre l'équation de la chaleur en 1D en fortran90 par différentes méthodes, et en premier lieu avec Euler explicite (j'ai également codé d'avance les subroutines RK2 et RK4). Les extrémités de la barre reste tout le temps à 0 degré et au temps t0, le reste de la barre est à 5 degrés. Par conséquent, la barre devrait refroidir, or le programme montre que la température augmente sans s'arrêter. J'ai beau chercher, je ne vois pas où est l'erreur (même en changeant les conditions initiales).
Merci d'avance!

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
module shemaeuler
  implicit none
 
contains
 
 subroutine fonct(Un,dt,dx,f,nx)
    implicit none
    real*8,dimension (nx),intent(in) :: Un
    real*8,dimension(nx),intent(out)::F
    real*8,intent(in)::dx,dt
    integer,intent(in)::nx
    integer::i
 
    F=0.0D0
    F(1)=0.0D0
    F(nx)=0.0D0
 
 
    do i=2,nx-1
       F(i)=(dt/dx**2)*(Un(i+1)-2*Un(i)+Un(i-1))+Un(i)
    end do
 
 
 
  ! print*,f                                                                                  
 
  end subroutine fonct
 
 
subroutine Euler( Un,Un1,dx,dt,nt,nx,F)
 
    implicit none
    real*8, intent(inout), dimension(nx) :: Un
    real*8, intent(out), dimension(nx) :: Un1
    real*8, intent(in) :: dt,dx
    real*8, dimension (nx),intent(inout):: F
    integer,intent(in)::nx,nt
    integer::i
 
 
 call fonct(Un,dt,dx,f,nx)
 
       Un1=Un+dt*F
 
 
 
 
  end subroutine Euler
 
end module shemaeuler
 
Program prog1
 
  use shemaeuler
 
  implicit none
 
  real*8, dimension (:),allocatable :: Un, Un1
  real*8, dimension (:),allocatable:: F
  real*8, dimension(:,:), allocatable :: mat
  real*8::dx,dt
  integer::i,n,tf,nx,nt
  real*8::D
 
 D=1D0
  dx=0.1D0
  dt=0.001D0
  tf=1.0D0
nx=floor(1./dx)
  nt=floor(tf/dt)
 
  allocate(Un1(nx),Un(nx),F(nx),mat(nx,nt))
 
  Un=5.0D0
  Un(1)=0.D0
  Un(nx)=0.0D0
  Un1=0.D0
 
 do i=1,nt
 
    call Euler( Un, Un1,dx,dt,nt,nx,f)
 
      mat(:,i) = Un1
     Un = Un1
      print*,'Un1=', Un1
 
  end do
end Program prog1

[xflr6]

Bonjour,

My guess:
Dans la fonction Euler vous résolvez l'ed temporelle Un1=Un+dt*F, mais vous avez déjà fait une partie du travail dans fonc avec :
F(i)=(dt/dx**2)*(Un(i+1)-2*Un(i)+Un(i-1))+Un(i).

Cordialement,
xflr6

[bengir]

Bonjour,
Je m'excuse je n'avais pas vu votre réponse. En effet je me suis rendu compte de la faute en ajoutant le terme un(i). Par contre est ce une faute de multiplier par dt dans ma fonction ? Car si je ne multiplie pas les résultats sont faux...
De plus j'ai une question supplémentaire. Dans mon projet on demande d'étudier le comportement de la solution en fonction du coefficient de diffusion D. Dans ma discretisation, il se retrouve à multiplier le terme 2*un(i)*dt/dx2. Problème, lorsque je dis qu'il vaut 1 (comme s'il n'y étais pas), j'obtiens des résultats complètement faux. À quoi est ce du ? Je peux vous envoyer le ccode si vous le souhaitez.
Merci d'avance

[xflr6]

Bonjour,

Comme pour l'ajout de un(i), il ne faut appliquer dt qu'une fois. La partie Un1=Un+dt*F, c'est la partie intégration temporelle, autant la garder dans la routine Euler, et enlever le un(i) et le dt dans la routine fonc.
Si j'ai bien compris et si vous gardez les deux dt, ce n'est plus l'équation de la chaleur qui est résolue.

Pour la seconde question, le coefficient de diffusion multiplie l'opérateur Laplacien, donc dans ce cas simple, tout ce qui est multiplié par 1/dx**2 doit l'être aussi par D, i.e. D/(dx*dx). Pas uniquement 2*un(i)*dt/dx2. Après, il y a aussi des conditions de stabilité pour D donné dans le choix des pas en temps et en espace (D*dt/dx**2 plus petit que 1 typiquement), mais je pense que le problème auquel vous êtes confronté est lié à la question précédente.
Un coup d'oeil sur le code pourrait être utile en effet.

Cordialement,
xflr6

[bengir]

Merci, ça marche maintenant.
Puis-je vous poser une autre question sur le même thème svp?
Je veux résoudre l'équation avec la méthode implicite. Donc il me faut résoudre Un=A*Un+1. A étant de grande dimension, j'utilise la factorisation LU pour résoudre ce problème. Là encore j'ai un souci. Ma factorisation est bonne (faite en TP et corrigée) mais mes résultats sont faux. Je pense que ça vient de ma matrice A (mon prof ne voit pas d'erreurs majeures mais soupçonne que ça vient de là). Quand je demande de l'afficher je me retrouve avec en 1ere ligne : 1 -alpha 0....0 et en seconde 0 1+2alpha -alpha. Or je devrais avoir :
1 0.....................0
-alpha 1+2alpha -alpha 0.....0
ce que je ne comprends pas c'est que quand je demande d'afficher seulement les 2 premières lignes elles sont justes! D'où est ce que ça peut venir?
Je vous joints le code.
De plus, je n'ai pas fait intervenir le coefficient D. Si je veux le mettre il suffit de multiplier A par D c'est exact? Car là encore les résultats sont faux.
Merci d'avance.

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
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
196
197
198
199
200
201
202
203
204
205
206
207
208
module fonctions
 
  !module principal contenant les subroutines
  implicit none
 
contains
 
!**************************************************************************
!Résolution d'un système linéaire triangulaire supérieur
 
subroutine Resol_Triang_Sup(n,A,B,X)
 
  implicit none
 
  integer, intent(in) :: n
  real*8, dimension(n,n),intent(in) :: A
  real*8, dimension(n), intent(in) :: B
  real*8, dimension(n), intent(out) :: X
  integer :: i, j
  real*8 :: somme
 
 
  !initialisation de X et somme
  X = 0
  somme = 0
 
  !algorithme de résolution
  X(n) = B(n)/A(n,n)
 
  Do i = n-1,1,-1
     somme = B(i)
     Do j = n, i, -1 
        somme = somme - A(i,j)*X(j)
     End Do
     X(i) = somme/A(i,i)
  End Do
 
end subroutine Resol_Triang_Sup
 
!*****************************************************************************
!Résolution d'un système linéaire triangulaire inférieur
 
subroutine Resol_Triang_Inf(n,A,B,X)
 
  implicit none
 
  integer, intent(in) :: n
  real*8, dimension(n,n),intent(in) :: A
  real*8, dimension(n), intent(in) :: B
  real*8, dimension(n), intent(out) :: X
  integer :: i, j
  real*8 :: somme
 
 
  !initialisation de X et somme
  X = 0
  somme = 0
 
  !algorithme de résolution
  X(1) = B(1)/A(1,1)
 
  Do i = 2, n
     somme = B(i)
     Do j = 1, i-1
        somme = somme - A(i,j)*X(j)
     End Do
     X(i) = somme/A(i,i)
  End Do
 
end subroutine Resol_Triang_Inf
!**************************************************************************
 
!Factorisation LU
 
subroutine LU(n,A,L,U)
 
  implicit none
 
  integer, intent(in) :: n
  real*8, dimension(n,n), intent(in) :: A
  real*8, dimension(n,n), intent(out) :: L, U
  integer :: i, j, k
 
  U = A
  L = 0
 
  Do i = 1, n
     L(i,i) = 1
  End Do
 
  Do i = 1, n
     Do j = i+1, n
        L(j,i) = U(j,i)/U(i,i)
        U(j,i) = 0
        Do k = i+1, n
           U(j,k) = U(j,k) - L(j,i)*U(i,k)
        End Do
     End Do
  End Do
 
end subroutine LU
!***************************************************************************
!Resolution LU
subroutine Resolution(n,L,U,B,X)
 
  implicit none
 
  integer, intent(in) :: n
  real*8, dimension(n,n), intent(in) :: L, U
  real*8, dimension(n), intent(in) :: B
  real*8, dimension(n), intent(out) :: X
  real*8, dimension(n) :: Y
  integer :: i, j
 
  call Resol_Triang_Inf(n,L,B,Y)
  call Resol_Triang_Sup(n,U,Y,X)
 
 
end subroutine Resolution
 
!***************************************************************
 
!ecriture dans un fichier
 
subroutine WriteGnuPlot(nom_fichier, v)
    implicit none
    character(len=*) :: nom_fichier
    real*8, dimension(:,:) :: v
    integer :: i,j
    open(unit = 11, file = nom_fichier, form = 'formatted', &
         status = 'unknown', action = 'write')
    do i = 1, size(v,1)
       do j=1,size(v,2)
          write(11, * )i, j, v(i,j)
       end do
       write(11, *)
    end do
    close(11)
  end subroutine WriteGnuPlot
 
end module fonctions
!****************************************************************************
 
 
!programme principal
program eulerimplicit
  use fonctions
  implicit none  
  real*8,dimension(:,:),allocatable::A,L,U,mat
  real*8,dimension(:),allocatable::Un,Un1,y
  real*8::Beta1,alpha,dt,dx,tf,D
  integer::i,j,nx,nt
  D=0.01
  tf=1.0D0
  dt=0.1d0
  dx=0.1d0
 
  nx=floor(1./dx)
  nt=floor(tf/dt)
 
  allocate(A(nx,nx))
  allocate(L(nx,nx))
  allocate(U(nx,nx))
  allocate(y(nx))
  allocate(Un(nx))
  allocate(Un1(nx))
  allocate(mat(nx,nt))
  Beta1=(1+2*(dt/dx**2))
  alpha=dt/dx**2
  A(1,1)=1.d0
  A(nx,nx)=1.d0
  Un(1)=0.d0
  Un(nx)=0.d0
 
  do i=2,nx-1
     A(i,i)=Beta1
     A(i,i+1)=-alpha
     A(i,i-1)=-alpha
  end do
 
 
  Un(1)=0.d0
  Un(nx)=0.d0
  do j=2,nx-1
     Un(j)=5.d0
  end do
  open(unit=24,file='resultatsimplicistp')
  print*,A
  do i=1,nt
     call LU(nt,D*A,L,U) !on multiplie la matrice A par le coefficient D'
     call Resolution(nt,L,U,Un,Un1)
     mat(:,i) = Un1
     Un=Un1
     ! print*,'la matrice L est :',L
     ! print*,'la matrice U est :',U
     print*, 'la matrice Un1 est :',Un1
  end do
  call writegnuplot('implicit.dat',mat)
  close(24)
 
  deallocate(A)
  deallocate(L)
  deallocate(U)
  deallocate(y)
  deallocate(Un)
  deallocate(Un1)
 
end program eulerimplicit

[bengir]

Finalement j'ai trouvé mon erreur. Elle venait de mes coefficients alpha et beta. Lorsque je mets un tf=1s, j'obtiens sensiblement les mêmes résultats qu'avec la methode explicite. Par contre si je mets tf=10. La compilation se fait mais lors de l'éxecution le message suivant s'affiche :
*** Error in `./test2': free(): invalid next size (normal): 0x0000000000f1cf10 ***
======= Backtrace: =========
/lib64/libc.so.6(+0x75a4f)[0x7f2d6988fa4f]
/lib64/libc.so.6(+0x7cd78)[0x7f2d69896d78]
./test2[0x4014aa]
./test2[0x4048d4]
./test2[0x40536b]
/lib64/libc.so.6(__libc_start_main+0xf5)[0x7f2d6983bd65]
./test2[0x400a59]
======= Memory map: ========
00400000-00407000 r-xp 00000000 00:28 3759416062 /net/malt/m/bgiralda/TER/Programmes/programmes2/test2
00606000-00607000 r--p 00006000 00:28 3759416062 /net/malt/m/bgiralda/TER/Programmes/programmes2/test2
00607000-00608000 rw-p 00007000 00:28 3759416062 /net/malt/m/bgiralda/TER/Programmes/programmes2/test2
00f14000-00f35000 rw-p 00000000 00:00 0 [heap]
7f2d6981a000-7f2d699ce000 r-xp 00000000 08:02 1442501 /usr/lib64/libc-2.18.so
7f2d699ce000-7f2d69bcd000 ---p 001b4000 08:02 1442501 /usr/lib64/libc-2.18.so
7f2d69bcd000-7f2d69bd1000 r--p 001b3000 08:02 1442501 /usr/lib64/libc-2.18.so
7f2d69bd1000-7f2d69bd3000 rw-p 001b7000 08:02 1442501 /usr/lib64/libc-2.18.so
7f2d69bd3000-7f2d69bd8000 rw-p 00000000 00:00 0
7f2d69bd8000-7f2d69c13000 r-xp 00000000 08:02 1443401 /usr/lib64/libquadmath.so.0.0.0
7f2d69c13000-7f2d69e12000 ---p 0003b000 08:02 1443401 /usr/lib64/libquadmath.so.0.0.0
7f2d69e12000-7f2d69e13000 r--p 0003a000 08:02 1443401 /usr/lib64/libquadmath.so.0.0.0
7f2d69e13000-7f2d69e14000 rw-p 0003b000 08:02 1443401 /usr/lib64/libquadmath.so.0.0.0
7f2d69e14000-7f2d69e29000 r-xp 00000000 08:02 1441807 /usr/lib64/libgcc_s-4.8.3-20140911.so.1
7f2d69e29000-7f2d6a028000 ---p 00015000 08:02 1441807 /usr/lib64/libgcc_s-4.8.3-20140911.so.1
7f2d6a028000-7f2d6a029000 r--p 00014000 08:02 1441807 /usr/lib64/libgcc_s-4.8.3-20140911.so.1
7f2d6a029000-7f2d6a02a000 rw-p 00015000 08:02 1441807 /usr/lib64/libgcc_s-4.8.3-20140911.so.1
7f2d6a02a000-7f2d6a12f000 r-xp 00000000 08:02 1442509 /usr/lib64/libm-2.18.so
7f2d6a12f000-7f2d6a32f000 ---p 00105000 08:02 1442509 /usr/lib64/libm-2.18.so
7f2d6a32f000-7f2d6a330000 r--p 00105000 08:02 1442509 /usr/lib64/libm-2.18.so
7f2d6a330000-7f2d6a331000 rw-p 00106000 08:02 1442509 /usr/lib64/libm-2.18.so
7f2d6a331000-7f2d6a450000 r-xp 00000000 08:02 1443404 /usr/lib64/libgfortran.so.3.0.0
7f2d6a450000-7f2d6a650000 ---p 0011f000 08:02 1443404 /usr/lib64/libgfortran.so.3.0.0
7f2d6a650000-7f2d6a651000 r--p 0011f000 08:02 1443404 /usr/lib64/libgfortran.so.3.0.0
7f2d6a651000-7f2d6a653000 rw-p 00120000 08:02 1443404 /usr/lib64/libgfortran.so.3.0.0
7f2d6a653000-7f2d6a673000 r-xp 00000000 08:02 1442918 /usr/lib64/ld-2.18.so
7f2d6a843000-7f2d6a847000 rw-p 00000000 00:00 0
7f2d6a870000-7f2d6a872000 rw-p 00000000 00:00 0
7f2d6a872000-7f2d6a873000 r--p 0001f000 08:02 1442918 /usr/lib64/ld-2.18.so
7f2d6a873000-7f2d6a874000 rw-p 00020000 08:02 1442918 /usr/lib64/ld-2.18.so
7f2d6a874000-7f2d6a875000 rw-p 00000000 00:00 0
7ffff37ad000-7ffff37ce000 rw-p 00000000 00:00 0 [stack]
7ffff37e3000-7ffff37e5000 r--p 00000000 00:00 0 [vvar]
7ffff37e5000-7ffff37e7000 r-xp 00000000 00:00 0 [vdso]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0 [vsyscall]

Program received signal SIGABRT: Process abort signal.

Backtrace for this error:
#0 0x7F2D6A34A497
#1 0x7F2D6A34AADE
#2 0x7F2D6984F8EF
#3 0x7F2D6984F877
#4 0x7F2D69850F67
#5 0x7F2D6988FA53
#6 0x7F2D69896D77
#7 0x4014A9 in __fonctions_MOD_resolution at neweulerimp.f90:112
#8 0x4048D3 in eulerimplicit at neweulerimp.f90:192
Abandon (core dumped)

Qu'est ce qu'il signifie?