Discussion :

[jimino]

Bonjour à tous, j'essaye de résoudre l'équation de la chaleur en dimension 1 sur un segment [a,b] par Fortran. J'ai réussi à faire un programme fonctionnel qui me donne à la fin un vecteur mais:
- J'obtiens à la fin encore le vecteur de condition initiale, pourquoi?
- Comment enregistrer à la fin mes valeurs de mon vecteur U? ( A la place je les copie colle du terminal puis je les trace dans un éditeur graphique )
- Comment utiliser un vecteur complexe en entrée? ( Par exemple, exp(ix) quelquechose ).

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
program  heatequation
implicit none
 
!Iteration variables
integer :: i
integer :: k
integer :: j
 
!Numerical Parameters
 
real, parameter :: xa=0.
real, parameter :: xb=15.
real, parameter :: dx =0.01
real, parameter :: dt=0.1
integer, parameter :: N=(xb-xa)/dx+1 
real, dimension(N,N) :: A
integer, parameter :: T=50
integer :: TT
real, dimension(N) :: U0 
real, dimension(N) :: U
real, dimension(N) :: V
real :: PI=3.1415926
 
!Initial condition 
do i=1,N
   U(i)=SIN((i-1)*PI*dx/xb)
end do
 
!Laplacian matrix
do j=1,N
do i=1,N
   A(i,j)=0
end do
end do
do k=1,N-1
   A(k,k)=-2
   A(k+1,k)=1
   A(k,k+1)=1
end do
A(N,N)=-2
 
!Solving PDE
TT=T/dt
do i=1,TT
   U=U+dt*matmul(A,U)
end do
 
print*,'U at time T is equal to',U
end program heatequation

[François L.]

Bonjour,

je vais commencer par répondre aux questions , et ferai ensuite quelques remarques:

* J'ai lancé le code, on n'obtient pas le vecteur initial après 50 unités de temps. Les valeurs ont changé. Ceci dit, quelles sont les CL ?
Parce que là, la première & la dernière ligne de A me paraissent suspectes

Par ailleurs, il manque une division par dx**2 du dt dans la formule de mise à jour de U.

* Il faut les écrire dans un fichier. Par exemple:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
OPEN(10,FILE="coincoin.res")
DO i=1,N
  WRITE(*,*) U(i)
ENDDO
CLOSE(10)

* Le fortran reconnaît le type COMPLEX

Quelques remarques:

* je préfère définir N d'abord & en déduire dx. Pour dt, je préfère une boucle du type

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
 DO WHILE t<=tfinal
   ...
  t=t+dt
ENDDO

* Quand je fais de la résolution d'EDP, je mets mes réels en double précision (ça dépend du compilateur, mais souvent on déclare REAL(8) au lieu de REAL )

* A peut être initialisée à 0 en écrivant simplement A=0.0d0

* Attention tout de même avec le schéma explicite, le dt/dx**2 me paraît très grand

Voila, hope this helps

[jimino]

Super, merci beaucoup! En effet, j'avais oublié le dx. Merci!
Matlab est-il compatible avec le fichier en sortie? (de sorte que je puisse le tracer graphiquement).

[François L.]

A priori oui, le programme de calcul écrira par défaut dans un fichier texte, qu'un load Matlab devrait pouvoir lire. Il y a peut-être la question des notations double précision du Fortran (i.e. 0.3D2 ), je ne me souviens plus si Matlab les digère. Mais de façon générale, j'utilisais souvent Matlab pour post-traiter mes résultats de code Fortran. Et je ne dois pas être le seul