Bonsoir,

J'ai réalisé un TP implémentant un schéma aux différences finies avec euler explicite, implicite et crank nicolson.
Le schéma est le suivant:
u(x,t)+u(x,t)+c=f(x,t)
u(0,t)=u(1,t)=0 pour tout t
u(x,0)=u0(x)

Voici le code pour Euler Explicite. Je dois pouvoir vérifier les ordres de convergences expérimentalement.
Je n'ai aucun souci pour vérifier l'ordre 2 en espace (pas divisé par 2, erreur par 4). Mais pour l'ordre 1 en temps, je n'y arrive pas.
Je me suis demandé si la condition de stabilité n'était pas un obstacle à cela.
Je dois utiliser des fonctions solutions non polynomiales en espace et en temps (ce qui est la cas ici).
Voilà le code, je suis bloqué, bloqué. Autre chose, si je ne respecte pas les conditions de stabilité mon erreur devrait exploser et non pas affiché l'erreur machine comme ici. Mon vecteur erreur est composé de nan, il y a donc un souci.

Merci d'avance pour votre aide

Voilà le code
Code:
Program eulerexp
Implicit None
Integer, parameter:: N=64, I=1000, c=1
Real, parameter :: T=2.0
REAL :: pasespace=1./I
Integer :: j, k
real, dimension (0:I,1:2) :: U
REAL, DIMENSION (0:I):: erreur

real :: pastemps
pastemps=T/N

Do j=0,I
U(j,1)=G(j*pasespace,0.0)
end do
! ordre 1 en temps et 2 en espace
! test espace: N=4000 I=10--I=20 Ok
!test temps: impossible à visualiser ? la stabilité est un obstacle?
! /!\ attention à la stabilité ici I*I<=N/4
!schéma numérique de euler explicite

do k=0,N-1
U(0,2)=0
U(I,2)=0

do j=1,I-1

U(j,2)=U(j,1)+ pastemps*(F(j*pasespace,k*pastemps)+((U(j+1,1)-2.0*U(j,1)+U(j-1,1))/(pasespace*pasespace))-c*((U(j+1,1)-&
U(j-1,1))/(2.0*pasespace)))

end do
U(:,1)=U(:,2)


end do

!calcul de l'erreur
do j=0,I
erreur(j)=U(j,1)-G(j*pasespace,T)
end do


print*, maxval(abs(erreur))


CONTAINS
	SUBROUTINE AFFICHE (A)
!
  REAL, DIMENSION (:,:) :: A
  INTEGER :: M,N,I,J
  CHARACTER*6 :: FMT='(F7.3)'
  M=SIZE(A,1)
  N=SIZE(A,2)

!
  DO I = 1,M
     DO J = 1,N-1
        WRITE (6,FMT,ADVANCE='NO') A(I,J)
     END DO
     WRITE (6,FMT,ADVANCE='YES') A(I,N)
  END DO
END SUBROUTINE AFFICHE

FUNCTION G(x,t)
Real :: x, G, t
G=exp(-t)*sin(2*3.14159*x)
End Function G

FUNCTION F(x,t)
Real ::x, F,t
F=-exp(-t)*sin(2.0*3.14159*x)+4.0*3.14159*3.14159*exp(-t)*sin(2.0*3.14159*x)+c*2.0*3.14159*exp(-t)*cos(2.0*3.14159*x)
END function F


end program eulerexp