Salve a tutti ! ho scritto una subroutine di modulo in fortran 90 che risolve la equazione a derivate parziali di burger , ora il mio dubbio e` questo ! la U e` calcolata in funzione di una variabile spaziale e una temporale !

Per testare il programma chiedo al programma di scrivermi su file il vettore X (relativo al dominio spaziale) e una colonna del vettore U da teoria il plot dovrebbe essere una semplice onda sinusoidale .. ma il problema e` che il programma calcola (e scrive su file) bene i valori della U(i , J=fissato) tranne gli ultimi 4 punti che tendono generalmente ad aumentare esponenzialmente !

Allego il codice (modulo e programma main )


codice:
  MODULE PARAM 
        IMPLICIT NONE 
        INTEGER, PARAMETER  :: SP = KIND(1.0)
        INTEGER, PARAMETER  :: DP = KIND(1.D0)
      END MODULE 
!-----------------------------------------------------------------------
      MODULE BURGER_2
       USE PARAM 
       IMPLICIT NONE
       
       REAL(DP), PARAMETER   :: x0=0. , xF= 5.
       REAL(DP), PARAMETER   :: t0=0. , tF= 1.0
       REAL(DP), PARAMETER   :: dx=0.01, dt=0.1
       INTEGER , PARAMETER   :: N= CEILING((xF-x0)/dx)
       INTEGER , PARAMETER   :: M= CEILING((tF-t0)/dt)
       REAL(DP), PARAMETER   :: nu = 2D-4 
       REAL(DP), PARAMETER   :: pi=3.14159265358979323846 




      CONTAINS 
!-----------------------------------------------------------------------
         SUBROUTINE BURGER_EXPLICIT (x,t,u) 
           USE PARAM 
           REAL(DP), DIMENSION(N+1)     :: x
           REAL(DP), DIMENSION(M+1)     :: t
           REAL(DP), DIMENSION(N+1,M+1) :: U
           REAL(DP)   :: r,k
           INTEGER    :: i,j,stat
            
          
          
           x(1:N+1)       = (/ ( x0+dx*(i-1) , i=1,N+1 ) /)       
           t(1:M+1)       = (/ (t0+ dt*(j-1) , j=1,M+1 ) /)
           U(1:N+1,1:M+1) = RESHAPE((/(((sin(2*pi*(i-1)/N)), i=1,N+1), j=1,M+1)/),(/N+1,M+1/)) 
             
           U(1,:)   = x0
           U(N+1,:) = xF


           U(:,1)   = t0
           U(:,M+1) = tF
           




           r = nu*DT/DX**2
           k = DT/(2*DX)
           
           OPEN(10,file='bur2.dat',STATUS='UNKNOWN',IOSTAT=stat)
           DO J=2,M
            DO I=2,N
              U(i,j+1) = r*u(i+1,j)+(1-2*r)*u(i,j)+r*u(i-1,j) - k*u(i,j)*(u(i+1,j)-u(i-1,j))
            END DO
           END DO 
              !IF(J .EQ. 4) THEN
              
              !WRITE(10,FMT='(2F20.5)') x(i), u(i,4)        !((u(i,j) , i=1,N+1), J=1,M+1)
              !END IF 


           !DO j=1,M+1  
           OPEN(10,file='bur1.dat',STATUS='UNKNOWN',IOSTAT=stat)
            IF (STAT .NE. 0) THEN 
             WRITE(6,FMT='(A)')'Error opening File ''bur1.dat'' (EXIT1)'
             STOP 
            ELSE     
              WRITE(10,FMT='(2F20.5)') (x(i), u(i,2), i=1,N)        !((u(i,j) , i=1,N+1), J=1,M+1)
              !WRITE(10,*)
              !WRITE(10,FMT='(2F20.5)') (x(i), u(i,2), i=1,N)        !((u(i,j) , i=1,N+1), J=1,M+1)
              WRITE(6,'(50("-"))') 
           !END DO 
            END IF


           RETURN   
         END SUBROUTINE 




!-----------------------------------------------------------------------
!-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
      END MODULE
e questo e` il main

codice:
      PROGRAM main        USE BURGER_2
       IMPLICIT NONE 


       REAL(DP) , DIMENSION (N+1) :: xx
       REAL(DP) , DIMENSION (M+1) :: tt
       REAL(DP) , DIMENSION (n+1,M+1) :: UU
      
       CALL burger_EXPLICIT(xx,tt,UU)


        
       STOP 
      END