Discussion :

[feynman]

Bonjour,
je cherche des programmes ecrit en fortran en difference finis pour les equations elliptiques ou paraboliques , est ce que quelqu un a des choses comme ca?
merci

[genteur slayer]

dans tous les cours d'analyse numérique....

enfin en tout cas tous ceux que j'ai lu donne des exemples en fortran....

[feynman]

Alors peut tu me donner un exemple qui se trouve sur le net? ou que tu posede?
merci

[genteur slayer]

http://goatin.univ-tln.fr/

c'est le site d'une de mes prof de maths, maintenant je lui fais ses TD (on comprend à voire sa photo pourquoi j'aime les maths maintenant) mais dans la rubrique teatching y a un cours d'ana num...

cependant pour les équation elliptique on utilise plutot un schéma eléments ou volumes finis car les diférences finies ne sont pas assez stable

[feynman]

merci genteur slayer ,
Moi je connais l analyse numerique, mais je cherche des programmes test en fortran pour la methode de difference finies , et non pas des cours theoriques en ana num
merci de nouveau

[genteur slayer]

bin c simple tu tape quasi texto la théorie dans un source fortran.... enfin c la base.... le cours de paola, les tp, on a fait que ça......

[feynman]

J ai essaye' sur les cours de Poala, y a pas des codes difference finie , mais volumes finies , moi je m interesse sur les difference finie,
merci pour ton aide

[genteur slayer]

dérivée, schémas avant de la fonction u:

u'~=(u_i+1 - u_i)/dl

schéma arrière:

u'~=(u_i - u_i-1)/dl

schéma centré

u'~=(u_i+1 - 2 u_i + u_i-1)/2dl

dérivée seconde:

u''~=(u_i+1 - 2 u_i + u_i-1)/dl^2

je vois toujours pas ce qu'il y a de compliqué là dedans

LEs diférences finies sont une manière de discrétiser tes équation donc à chaque équation corespond un 'programme' de résolution.....
c'est tellement tout bête que voilà....

imaginons que tu doivent résoudre l'équation de la chaleur super simplifié, l'équation donne en 2d: d2T/dx2+d2T/dy2=0 on se place sur une plaque carré, avec deux coté adjacent à 20° et les deux autres à 50°. on discrétise facilement en diférences finies, on dit que T(x,y)~=T(xi,yj)=Tij (note: on peut aussi numéroter avec un seul indice)on va dire que l'on a N points par coté,et h=1/N on a donc: dx=dy=h
(ti+1,j + ti,j+1 -4 tij +ti,j-1 +ti-1,j)/h² = 0

comme c bizare, tu transforme l'écritue en AX=b systeme matricielle où la matrice A est une matrice pentadiagonalegenre:

A= 1/h²[-4 1 0 ....0 1 0 ...........]
|1 -4 1 0 ... 0 1 0 ........
|0 1 4
|.
|.
|1
et X est un vecteur où l'on a Xi+N*j=Tij

et dans le b, c'est là que l'on calle nos conditions aux limites (je détaille pas c chiant)

tu résout un systèeme linéaire pas la méthode que tu veux (complètement in dépendant du schéma de diférence finie)

ET voili voilou

[phy4me]

bonjour

Je peux vous passer des codes que j'ai ecris pour les équations hyperbolique , et c'est facile de les utilisés pour votre cas, ils suffit juste de bien écrire les équations discétisés

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
! Résolution numérique de l'equation de transport en utilisant la méthode  !!
! de differences finis et le schéma UPWIND du 1er ordre                    !!
! Auteur Elmiloud CHAABELASRI                                              !!  
! Chaabelasri@gmail.com                                                    !! 
!===========================================================================     
      program TRANS_UPWIND
!==== Declaration des parametres============================================
      real  nx,nt,sigma,dx
 
      parameter(nx=100.0) !Nombre de pas dans espace 
 
	dimension up1(nx+1),up2(nx+1),x(nx+1),ex(nx+1)
      dx=0.02             !Pas d'espace 
      sigma=1.d0           !CFL
      nt=1.0/(sigma*dx)   !Nombre de pas dans le temps
	dt=sigma*dx
 
 
	write(*,*) 'nx=',nx
	write(*,*) 'nt=',nt
	write(*,*) 'dt=',dt
!==== Conditions initiales==================================================
!    # tranport d'une discontinuité
      up1(1)=1.0
	do i=2,nx+1
	up1(i)=0.0
	enddo
 
!==== Schema UPWIND du 1er ordre============================================
      do 100 t=1,nt
!----->
	 do i=2,nx
	  up2(i)=(1-sigma)*up1(i)+sigma*up1(i-1)
	 enddo
 
       up2(1)  =1.0
	 up2(101)=0.0
 
 
       do i=1,nx+1
	  up1(i)=up2(i)
	 enddo
 
!----->
100   continue
!==== Enregistrement des resultats==========================================
      open(64,FILE='sol sumer CFL 1.0')
	open(65,FILE='solution exact')
 
***** solution exacte
 
	do i=1,nx/2.d0
	ex(i)=1.d0
	enddo
 
	do i=(nx/2.d0)+1,nx+1
	ex(i)=0.0
	enddo
 
******
 
 
      x(1)=0.
	do i=2,nx+1
	x(i)=x(i-1)+dx
	enddo
 
	do i=1,nx+1
	write(64,999)x(i),up1(i)
	write(*,999) x(i),up1(i)
	end do
 
	do i=1,nx+1
	write(65,999)x(i),ex(i)
	end do
 
999	format(2f12.3)
 
      end

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
c burgers + Lax-Frederichs
c 
c Chaabelasri@gmail.com
 
       program burrglaxfreid
 
       dimension u1(1001),u(1001)   
 
	 xmax=2.d0
	 xmin=0.d0
	 dx=0.02 
 
	 cfl=0.1 
	 umax=2.d0
	 t=1.0   
 
	   nx=(xmax-xmin)/dx
 
	   dt=(cfl*dx)/umax
	   nt=t/dt
 
	   tau=dt/dx
 
       write(*,*) nt
	 write(*,*) nx
 
c Solution initiale
       u(1)=2.0
	 do i=1,nx+1
	 u(i)=1.0
	 enddo
 
c***** Boucle temps *****************
      do j=1,nt
	   do i=2,nx+1
 
         ! 1er ecriture
          u1(i)=(0.5)*(u(i+1)+u(i-1))-(dt/(4*dx))*(u(i+1)**2-u(i-1)**2)
 
	   ! 2eme ecritute	    
	   !fd=(1.0/4.0)*(u(i)**2+u(i+1)**2)-(dx/(2*dt))*(u(i+1)-u(i  ))
	   !fg=(1.0/4.0)*(u(i)**2+u(i-1)**2)-(dx/(2*dt))*(u(i  )-u(i-1))
	   !u1(i)=u(i)-(dt/dx)*(fd-fg)
 
 
 
	   enddo
 
c Condition aux limites
         u1(1)   =2.0
	   u1(nx+1)=1.0
 
C Initialisation
         do i=1,nx+1
	   u(i)=u1(i)
	   enddo
 
      enddo
c***** Fermeture boucle temps ******
 
c Lecture des solutions
 
	open(64,file='solution calculée cfl=0.1')
 
	xi=-dx     
	do i=1,nx+1
	xi=xi+dx
	write(64,999)xi,u1(i)
	write(*, 999)xi,u1(i)
	end do
 
 
999	format(2f12.4)
	end

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
!
!  Equation de convection lineaire avec terme source constant par schema de LF
!  
!  Chaabelasri@gmail.com
!
!
      Program ECLTCLF
 
      implicit double precision(a-h,o-z)  
	dimension u1(101),u2(101),z (101),x (101)
 
	real nx,t,dx,CFL,dt,nt
 
	nx =100
	t =2.0 
	dx =0.2
	CFL=1.0
 
	dt=CFL*dx 
	nt=t/dt
 
	write(*,*) 'nt=',nt,'-- dt=',dt,'-- CFL=',CFL    
 
!--------------------------------------
! Condition initials
!--------------------------------------      
 
       do i=1,50
	  u1(i)=6.0
	 enddo
 
	 do i=51,101
	  u1(i)=1.0
	 enddo
!----
 
       do i=1,50
	  z(i)=0.0
	 enddo
 
	 do i=51,101
	  z(i)=0.1
	 enddo
 
!--------------------------------------
! Boucle de temps
!--------------------------------------
 
      do 100 j=1,nt
 
!---> Schema Lax-Friedrichs
      do i=2,nx
	 r=(0.5)*(u1(i+1)+u1(i-1))
	 u2(i)=r-(dt/(2*dx))*(u1(i+1)-u1(i-1)+z(i+1)-z(i-1))
      enddo
 
!---> Conditions aux limite
       u2(1)  =6.0
	 u2(nx+1)=1.0
 
!---> Changement des matrices
      do i=1,nx+1
	u1(i)=u2(i)
	enddo
!---> Remonter vers le début de la boucle de temps
100   Continue
 
!--------------------------------------
! Lecture des résultats
!--------------------------------------
 
      open(UNIT=65,FILE='CFL=1.0',STATUS='new')
 
      x(1)=0.
	do i=2,nx+1
	x(i)=x(i-1)+dx
	enddo
 
	do i=1,nx+1
	write(65,999)x(i),u1(i)
	write(*,999) x(i),u1(i)
	end do
 
 
999	format(2f12.4)	
	end

J'espère que ça vous sera utile

conseil : essaye d'écrire ton propre code, biensur ces code vous permet juste de voir les technique de programmation, mais la bonne maitrise de programmation se developpe facilement par la manipulation --> bonne courage

[genteur slayer]

ça ressemble bcp à des TP d'école.... qu'est-ce que j'ai pu me faire chié avec l'equation de la chaleur ou celle de burgers!!!!

[feynman]

Merci phy4me pour votre aide, je voit que tu traite les probleme unidimentionel , as tu des exemples sur deux dimensions?
Moi je connais les methodes numeriques, mais je cherche les codes seulement

[phy4me]

d'acord, je vais vous passer un code écris par un professeur connu dans le domaine de la simulation numérique, c'est le professeur Toro, c'est un code en 2D (maillage rectangulaire sturcuté) qui résout les équations des écoulements d'eau peu profonde - shallow water - ( les équations de Saint venant ).
pour l'éxecution vous devez besoin d'un fichier de donnée je vais essayer de l'attaché avec cette réponse

bon courage

NB : ce code utilise la méthode des volumes finis

[feynman]

Ma demande est claire : exemple de code fortran pour la methode difference finies en 2D

[phy4me]

Alors prends tes équations discétisées, et allez à la soubroutine qui calcule les états ( les inconnus ) et essaye de la modiffiée, je pense qu'il y'a pas plus simple que ça tant que vous maitrisez l'analyse numérique comme vous l'avez mentioné.
vous n'allez jamais trouvez éxactement ce que vous chercher, un petit effort est demandé,

[feynman]

Alors quelles sont ces subroutines? ou je peux les trouver?

[phy4me]

La soubroutine principale de calcul est : SWEEPS ( appelée par le PP)
dans cette soubroutine il ya une appelation d'autres soubrutines, mais le volumes finis apparait dans la soubroutine RIEMAN ou les états gauche et droite sont calculés approximativement par un solveur Riemann approché SLLH, c'est la ou le principale changement devait etre,

[feynman]

Puisque je suis debutant, comment je les utilises ces subroutines? comment et ou je les fait un call?
merci

[FR119492]

Trois remarques:

Primo: Comme un autre membre te l'a déjà signalé, dans la plupart des cas, la méthode des différences finies est avantageusement remplacée par celle des éléments finis, en particulier lorsque la frontière du domaine n'est pas constituée exclusivement de segments de droites parallèles aux axes de coordonnées.

Secundo: Puisque tu connais la méthode pour construire le système linéaire, programme-la. Ensuite, tu vas sur www.netlib.org et tu télécharges la bibliothèque LinPack. Pour résoudre ton système, tu utilises les routines DPBFA ou DPBCO, puis DPBSL. Les lignes de commentaires au début de chaque routine t'expliquent comment l'utiliser.

Tertio: Il manque une information essentielle dans ta question: combien de noeuds envisages-tu de définir: 50 ou 50'000? Du point de vue mathématique, ça ne change pas grand chose, mais du point de vue programmation, ça change tout, selon que tu dois travailler sur disque ou non.

Bonne chance
Jean-Marc Blanc

[feynman]

Merci Jean Marc pour tes explications, ca marche un peu mieux maintenant,
Je suis oblige' d utiliser la methode de difference finie vu que tout mon travail est sur le bord (pour plus de precision, moi je calcul le flux de la solution de l equation div(\gamma\nabla u)=0 et u=f sur le bord ). Alors la methode des elements finis ne sert a rien dans mon cas , sinon il y a freefem pour programmer la methode des element finie et c est plus facile que fortran.
Mon domaine c est le cercle, et c est facile de faire les difference finie la dessous. Mon seul probleme est que je ne connais pas fortran de tout, jusqu a le mois passe', j ai commence'a programmer en utilisant le fortran, donc je suis vraiment debutant

[FR119492]

Salut.
Si ton domaine d'intégration est un cercle, comme j'ai cru le comprendre en lisant ton dernier message, et si la répartition de gamma est de révolution, pourquoi n'écris-tu pas ton équation aux dérivées partielles en coordonnées polaires. Dans ce cas, les variables se séparent et tu obtiens une équation différentielle ordinaire. Ton problème devient alors beaucoup plus simple.
Bonne chance
Jean-Marc Blanc