!****************************************************************************!
!Programme de calcul de l'exercice 1 de la la feuille 6                      !
!Verification de l'ordre des schemas suivants                                !
!*schema d'Euler explicite                                                   !
!*schema de Taylor ordre 2                                                   !
!*schema de Taylor ordre 3                                                   !
!*schema de Runge-Kutta d'ordre 2                                            !
!****************************************************************************!

!REMARQUE: il vaut mieux compiler en double precision

program ordre
  implicit none
  real                    ::h           !le pas de temps
  real                    ::y1          !la solution exacte
  real                    ::y2,y3,y4,y5 !la valeur des differents schemas
  real                    ::e2,e3,e4,e5 !les erreurs des differents schemas
  real                    ::t           !le temps
  integer                 ::n           !nombre de pas de temps
  integer                 ::i           !l'increment de temps
  integer                 ::k           !l'increment de subdivision




  open(unit=10,file='ordre.dat',status='unknown')
  write(10,*) '#log(h)   euler  taylor2 taylor3 RK2'
  
  do k=4,14
     n=2**k
     
     !calcul du pas de temps
     h=4.*atan(1.)/n
     
     
     !initialisation
     y1=1.
     y2=1.
     y3=1.
     y4=1.
     y5=1.
     t=0.
     e2=0.
     e3=0.
     e4=0.
     e5=0.
     
     
     do i=1,n
        y2=y2+h*phi2(t,y2,h)
        y3=y3+h*phi3(t,y3,h)
        y4=y4+h*phi4(t,y4,h)
        y5=y5+h*phi5(t,y5,h)
        t=t+h
        y1=exp(1.-cos(t))
        e2=max(e2,abs(y2-y1))
        e3=max(e3,abs(y3-y1))
        e4=max(e4,abs(y4-y1))
        e5=max(e5,abs(y5-y1))
     enddo
     
     write(10,200) log(h),log(e2),log(e3),log(e4),log(e5)
  end do
  close(10)
  print*,'fin du programme'
  stop


200   format (5(E16.7,1x))


contains

  function phi2(t,y,h) result(phi)
    !schema d'Euler explicite
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)
    
  end function phi2


  function phi3(t,y,h) result(phi)
    !schema de Taylor, ordre 2
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)+0.5*h*(y*cos(t) + y*sin(t)*sin(t))
    
  end function phi3

  function phi4(t,y,h) result(phi)
    !schema de Taylor, ordre 3
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=y*sin(t)+0.5*h*(y*cos(t) + y*sin(t)*sin(t)) & 
         &+(h**2)/6.*(-y*sin(t) +3*y*sin(t)*cos(t)+y*(sin(t))**3)
    
  end function phi4
  

  function phi5(t,y,h) result(phi)
    !schema de Runge Kutta, ordre 2
    implicit none
    real            ::t,y,h
    real            ::phi
    
    phi=sin(t+0.5*h)*(y+0.5*h*sin(t)*y)
    
  end function phi5
end program ordre
